Matlab中二阶常微分方程的求解系统

Question

简介

我正在使用 Matlab 通过使用 ODE45 数值求解二阶常微分方程的系统来模拟一些动态系统。我从 Mathworks 中找到了一个很棒的教程（最后的教程链接）关于如何做到这一点。

在本教程中，方程组在 x 和 y 中是显式的，如下所示：

x''=-D(y) * x' * sqrt(x'^2 + y'^2)
y''=-D(y) * y' * sqrt(x'^2 + y'^2) + g(y)

以上两个方程的形式为 y'' = f(x, x', y, y')

---

问题

但是，我遇到了无法明确求解变量的方程组，如示例中所示。例如，其中一个系统具有以下一组 3 个二阶常微分方程：

y 双素数方程

y'' - .5*L*(x''*sin(x) + x'^2*cos(x) + (k/m)*y - g = 0

x 双素数方程

.33*L^2*x'' - .5*L*y''sin(x) - .33*L^2*C*cos(x) + .5*g*L*sin(x) = 0

一个素数是一阶导数 双素数是二阶导数 L、g、m、k 和 C 是给定的参数。

如何使用 Matlab 数值求解一组二阶常微分方程，其中二阶无法显式求解？

谢谢！

Answer 1

你的第二个系统有表格

a11*x'' + a12*y'' = f1(x,y,x',y')
a21*x'' + a22*y'' = f2(x,y,x',y')

您可以将其作为线性系统求解

[x'', y''] = A\f

或者在这种情况下明确使用克莱默规则

x'' = ( a22*f1 - a12*f2 ) / (a11*a22 - a12*a21)

y'' 相应地。

我强烈建议在代码中保留中间变量，以减少输入错误的机会并避免对相同表达式进行多次计算。

代码可能看起来像这样（未经测试）

function dz = odefunc(t,z)
    x=z(1); dx=z(2); y=z(3); dy=z(4);
    A = [  [-.5*L*sin(x),  1] ;  [.33*L^2,  -0.5*L*sin(x)]  ]
    b = [ [dx^2*cos(x) + (k/m)*y-g]; [-.33*L^2*C*cos(x) + .5*g*L*sin(x)] ]

    d2 = A\b

    dz = [ dx, d2(1), dy, d2(2) ]
end

Answer 2

是的，您的方法是正确的！ 我在下面发布以下代码：

%Rotating Pendulum Sym Main
clc
clear all;

%Define parameters

global M K L g C;
M = 1;
K = 25.6;
L = 1;
C = 1;
g = 9.8;


% define initial values for theta, thetad, del, deld
e_0 = 1;
ed_0 = 0;
theta_0 = 0;
thetad_0 = .5;
initialValues = [e_0, ed_0, theta_0, thetad_0];

% Set a timespan
t_initial = 0;
t_final = 36;
dt = .01;
N = (t_final - t_initial)/dt;
timeSpan = linspace(t_final, t_initial, N);

% Run ode45 to get z (theta, thetad, del, deld)
[t, z] = ode45(@RotSpngHndl, timeSpan, initialValues);

%initialize variables
e = zeros(N,1);
ed = zeros(N,1);
theta = zeros(N,1);
thetad = zeros(N,1);
T = zeros(N,1);
V = zeros(N,1);
x = zeros(N,1);
y = zeros(N,1);

for i = 1:N
    e(i) = z(i, 1);
    ed(i) = z(i, 2);
    theta(i) = z(i, 3);
    thetad(i) = z(i, 4);
    T(i) = .5*M*(ed(i)^2 + (1/3)*L^2*C*sin(theta(i)) + (1/3)*L^2*thetad(i)^2 - L*ed(i)*thetad(i)*sin(theta(i)));
    V(i) = -M*g*(e(i) + .5*L*cos(theta(i)));
    E(i) = T(i) + V(i);
end

figure(1)
plot(t, T,'r');
hold on;
plot(t, V,'b');
plot(t,E,'y');
title('Energy');
xlabel('time(sec)');
legend('Kinetic Energy', 'Potential Energy', 'Total Energy');

这里是 ode45 的函数句柄文件：

function dz = RotSpngHndl(~, z)
% Define Global Parameters
global M K L g C

A = [1, -.5*L*sin(z(3));
    -.5*L*sin(z(3)), (1/3)*L^2];

b = [.5*L*z(4)^2*cos(z(3)) - (K/M)*z(1) + g;
    (1/3)*L^2*C*cos(z(3)) + .5*g*L*sin(z(3))];

X = A\b;

% return column vector [ed; edd; ed; edd]
dz = [z(2);
    X(1);
    z(4);
    X(2)];