大神求解动力学微分方程,用matlab中ode45编程实现!

大神求解动力学微分方程,用matlab中ode45编程实现!

其中的参数可以随便设置,

参考代码：
 
% 参数定义
m1 = 1; m2 = 2; m3 = 3;
k1 = 4; k2 = 5; k3 = 6; k4 = 7;
P0 = 8; w = 9;
% 微分方程
A = [k1+k2 -k1 0; -k2 k2+k3 -k3; 0 -k3 k3+k4];
dx = @(t,x)[x(4:6); (-A*x(1:3)+[P0*sin(w*t);0;0])./[m1;m2;m3]];
% 初始条件（x4~x6 为 x1'~x3'）
x0 = zeros(6,1);
% 直接绘图
ode45(dx,[0 5],x0)
lstr = arrayfun(@(i){sprintf('x%i',i)},1:6);
legend(lstr{:})
legend(lstr{:},'Location','S','Orientation','horizontal')
 

再问： 在微分方程中加入一个阻尼矩阵怎么编程，还有程序里的@是什么意思，主要是matlab都不会。。。。
再答： 程序中的@ 表示匿名函数（anonymous function）。关于匿名函数的更多信息，请自行搜索相关介绍。 要想加入阻尼，除了定义常数C1-C3，另外需要修改一下匿名函数dx： % 参数定义
m1 = 1; m2 = 2; m3 = 3;
k1 = 4; k2 = 5; k3 = 6; k4 = 7;
P0 = 8; w = 9;
C1 = 1; C2 = 2; C3 =3;% 微分方程
A = [k1+k2 -k1 0; -k2 k2+k3 -k3; 0 -k3 k3+k4];
dx = @(t,x)[x(4:6); (-A*x(1:3)-x(4:6)./[C1;C2;C3]+[P0*sin(w*t);0;0])./[m1;m2;m3]];% 初始条件（x4~x6 为 x1'~x3'）
x0 = zeros(6,1);% 求解，直接绘图
ode45(dx,[0 5],x0)
lstr = arrayfun(@(i){sprintf('x%i',i)},1:6);
legend(lstr{:})
legend(lstr{:},'Location','S','Orientation','horizontal') 
再问： 如果阻尼不是对角阵，怎么编写匿名函数呢？如阻尼矩阵是
再答： 上次写错了，匿名函数   dx = @(t,x)[x(4:6); (-A*x(1:3)-x(4:6)./[C1;C2;C3]+[P0*sin(w*t);0;0])./[m1;m2;m3]];应为    dx = @(t,x)[x(4:6); (-A*x(1:3)-x(4:6).*[C1;C2;C3]+[P0*sin(w*t);0;0])./[m1;m2;m3]];如果C不是对角阵，就改成    dx = @(t,x)[x(4:6); (-A*x(1:3)-C*x(4:6)+[P0*sin(w*t);0;0])./[m1;m2;m3]]; 不举例子了，你自己试试吧。