Question

二元二阶微分方程组求数值解的matlab程序

Answer 1

方程或者初始条件有问题，请认真核实。

原因：如果按照现在的初始条件，则第一个微分方程各项为0，方程不起作用。而且就算法而言，二阶（最高阶）导数的系数为0，会导致数值积分算法无法启动。

 

参考代码大致如下：

m=1;
g=9.8;
k=100;
l0=1.1;
x0=[0.1 0 0 0];

% 定义x1=l, x2=l', x3=a, x4=a'
dx=#(t,x)[x(2); (m*x(1)*x(4)^2-m*g*cos(x(3))-k*(x(1)-l0))/m; ...
    x(4); (-2*m*x(1)*x(2)*x(4)+m*g*x(1)*sin(x(3)))/(m*x(1)^2)];
[t,x]=ode15s(dx,[0 10],x0);
lstr = { '\itl', '{\itl}''', '\alpha', '\alpha''' };
for i=1:length(lstr)
    subplot(2,2,i)
    plot(t, x(:,i));
    xlabel('Time')
    ylabel( lstr{i} )
end

注意：其中的#应为电子邮件的AT符，在MATLAB中表示匿名函数。因出现这个字符就会被识别为邮件地址，很可能回答无法正常显示（恶心的度娘，鄙视），所以，慎重起见，这里以#代替。

Answer 2

方程或者初始条件有问题，请认真核实。
原因：如果按照现在的初始条件，则第一个微分方程各项为0，方程不起作用。而且就算法而言，二阶（最高阶）导数的系数为0，会导致数值积分算法无法启动。

参考代码大致如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

m=1;
g=9.8;
k=100;
l0=1.1;
x0=[0.1 0 0 0];

% 定义x1=l, x2=l', x3=a, x4=a'
dx=#(t,x)[x(2); (m*x(1)*x(4)^2-m*g*cos(x(3))-k*(x(1)-l0))/m; ...
x(4); (-2*m*x(1)*x(2)*x(4)+m*g*x(1)*sin(x(3)))/(m*x(1)^2)];
[t,x]=ode15s(dx,[0 10],x0);
lstr = { '\itl', '{\itl}''', '\alpha', '\alpha''' };
for i=1:length(lstr)
subplot(2,2,i)
plot(t, x(:,i));
xlabel('Time')
ylabel( lstr{i} )
end

Answer 3

function dz=yueh(t,z) %保存为 yueh.m 文件
dz=[z(2);
10-z(1)-0.001*z(3);
z(4);
-0.0001*z(1)-0.001*z(3)];
-------------------------------------------------------
运行一下：
[t,z]=ode45('yueh',[0 1],[0 0 0 0]);
% 感觉初始值不完整啊,这里假设x(0)=θ(0)=x'(0)=θ'(0)=0
plot(z(:,1),z(:,3)) % 这里是 x-θ 图像

Answer 4

微分方程初值问题，用ode45可以解