Question

大神求解动力学微分方程，用matlab中ode45编程实现！谢谢

其中的参数可以随便设置，谢谢！

Answer 1

参考代码：

 

% 参数定义
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]];

 

不举例子了，你自己试试吧。