10自由度二阶微分方程组的MATLAB求解问题

fineshang

我正在分析一个振动问题，列出的是10自由度二阶微分方程组，形如：
[m]{X''}+[C]{X'}+[K]{X}={P}
其中[m]=diag{J1,J2,J3,J4,J5,J6,J7,J8,m,m}
    {X}=（φ1，φ2，φ3，φ4，φ5，φ6，φ7，φ8，X,Y）
    [C]为10*10阶的常数矩阵
    [K]为10*10阶的时变矩阵，其中每个元素均是时间t的函数，Kij=Kij(t)
    {P}=diag{Tin,0,0,0,0,0,0,0,0,-Tout}
想用数值解法（四阶龙格库塔法），借助matlab进行求解，最后输出
   （1）X与Y的轨迹图
   （2）Xsp(t)=Rs×φ2-Rp×φ3-Rc×φ7 与 X'sp(t)=Rs×φ'2-Rp×φ'3-Rc
×φ'7
   （3）Xcl(t)=Rc×(φ7-φ8) 与 X'cl(t)=Rc×(φ'7-φ'8)
   （4）Xsp(t) 与 X'sp(t)的相平面图
可是，我不会用matlab，请高手赐教，怎么编程求解，非常感谢！！

[ 本帖最后由 eight 于 2007-8-20 20:35 编辑 ]

咕噜噜

用最简单的ode45就可以
很多书上都有介绍

花如月

问题不小呀，等待做过这个问题的高手帖程序上来吧。不过希望渺茫哦，建议你一边等，一边自己做。不会的就学，要是都会了，就用不着学习了:lol

octopussheng

天哪，十自由度、二阶微分方程组，兄弟，你比我还狠，这个太复杂了，解方程不好解啊！

你要是用ode方法的话，就参考本板块的一个关于微分方程的精华贴，照着做就可以了，不过解这样的方程可是需要时间的啊！要耐心哦！

xjzuo

将原问题贴出来看看.
可以先试试ode45之类的函数．

crystallemon

下面的程序是9自由度，二阶微分方程的rk解法：第一部分是定义函数：
定义函数
function ydot=f(t,Y,P)
m1=700
m2=150;%
m3=150;%对
m4=1700;%
m5=800;%
m6=600;%
i1=4;%
i2=3;%
i3=3;%

r1=0.225;%
r2=0.28;%
r3=0.28;%
k1=7000000;
k3=600000;
k4=600000;
kp=2000000;
ee=80000000000;
aa=0.000133;
nline=2;
lline=20;
kr1=nline*ee*aa/lline;
kr2=kr1;
khead=800000;
kr1s=kr1*khead/(kr1+khead);
kr2s=kr2*khead/(kr2+khead);
km11=k1+kr1+kr2;
km22=kr1+kr1s+k3;
km33=kr2+kr2s+k4;
km44=k4;
km55=k3+kp;
km66=kp;
km77=(kr1+kr2)*r1*r1*1.01;
km88=(kr1+kr1s)*r2*r2;
km99=(kr2+kr2s)*r3*r3;
km17=(kr1-kr2)*r1;
km27=-kr1*r1;
km37=kr2*r1;
km18=kr1*r2;
km28=(-kr1+kr1s)*r2;
km78=kr1*r1*r2;
km19=-kr2*r3;
km39=(kr2-kr2s)*r3;
km79=kr2*r1*r3;
M=[m1 0 0 0 0 0 0 0 0;0 m2 0 0 0 0 0 0 0;0 0 m3 0 0 0 0 0 0;0 0 0 m4 0 0 0 0 0;0 0 0 0 m5 0 0 0 0;0 0 0 0 0 m6 0 0 0;0 0 0 0 0 0 i1 0 0;0 0 0 0 0 0 0 i2 0;0 0 0 0 0 0 0 0 i3];
K=[km11 -kr1 -kr2 0 0 0 km17 km18 km19;-kr1 km22 0 0 -k3 0 km27 km28 0;-kr2 0 km33 -k4 0 0 km37 0 km39;0 0 -k4 km44 0 0 0 0 0;0 -k3 0 0 km55 -kp 0 0 0;0 0 0 0 -kp kp 0 0 0;km17 km27 km37 0 0 0 km77 km78 km79;km18 km28 0 0 0 0 km78 km88 0;km19 0 km39 0 0 0 km79 0 km99];
C=zeros(9,9);
%C=0.01*[km11 -kr1 -kr2 0 0 0 km17 km18 km19;-kr1 km22 0 0 -k3 0 km27 km28 0;-kr2 0 km33 -k4 0 0 km37 0 km39;0 0 -k4 km44 0 0 0 0 0;0 -k3 0 0 km55 -kp 0 0 0;0 0 0 0 -kp kp 0 0 0;km17 km27 km37 0 0 0 km77 km78 km79;km18 km28 0 0 0 0 km78 km88 0;km19 0 km39 0 0 0 km79 0 km99];
I=eye(9,9);
A=[zeros(9,9),I;-inv(M)*K,-inv(M)*C];
ydot=A*Y+P;

crystallemon

第二部分：求解
clear;
clc;
t0=0;
Y0=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
h=0.01;
tN=5;
w=150;
P0=5000;
PA=100;
j=1;
t=t0:h:tN;
N=length(t);

k1=7000000;
k3=600000;
k4=600000;
kp=2000000;
ee=80000000000;
aa=0.000133;
nline=2;
lline=20;
kr1=nline*ee*aa/lline;
kr2=kr1;
khead=800000;
kr1s=kr1*khead/(kr1+khead);
kr2s=kr2*khead/(kr2+khead);
km11=k1+kr1+kr2;
km22=kr1+kr1s+k3;
km33=kr2+kr2s+k4;
km44=k4;
km55=k3+kp;
km66=kp;
km77=(kr1+kr2)*r1*r1*1.01;
km88=(kr1+kr1s)*r2*r2;
km99=(kr2+kr2s)*r3*r3;
km17=(kr1-kr2)*r1;
km27=-kr1*r1;
km37=kr2*r1;
km18=kr1*r2;
km28=(-kr1+kr1s)*r2;
km78=kr1*r1*r2;
km19=-kr2*r3;
km39=(kr2-kr2s)*r3;
km79=kr2*r1*r3;
M=[m1 0 0 0 0 0 0 0 0;0 m2 0 0 0 0 0 0 0;0 0 m3 0 0 0 0 0 0;0 0 0 m4 0 0 0 0 0;0 0 0 0 m5 0 0 0 0;0 0 0 0 0 m6 0 0 0;0 0 0 0 0 0 i1 0 0;0 0 0 0 0 0 0 i2 0;0 0 0 0 0 0 0 0 i3];
K=[km11 -kr1 -kr2 0 0 0 km17 km18 km19;-kr1 km22 0 0 -k3 0 km27 km28 0;-kr2 0 km33 -k4 0 0 km37 0 km39;0 0 -k4 km44 0 0 0 0 0;0 -k3 0 0 km55 -kp 0 0 0;0 0 0 0 -kp kp 0 0 0;km17 km27 km37 0 0 0 km77 km78 km79;km18 km28 0 0 0 0 km78 km88 0;km19 0 km39 0 0 0 km79 0 km99];
%C=zeros(9,9);
C=0.01*[km11 -kr1 -kr2 0 0 0 km17 km18 km19;-kr1 km22 0 0 -k3 0 km27 km28 0;-kr2 0 km33 -k4 0 0 km37 0 km39;0 0 -k4 km44 0 0 0 0 0;0 -k3 0 0 km55 -kp 0 0 0;0 0 0 0 -kp kp 0 0 0;km17 km27 km37 0 0 0 km77 km78 km79;km18 km28 0 0 0 0 km78 km88 0;km19 0 km39 0 0 0 km79 0 km99];
I=eye(9,9);
A=[zeros(9,9),I;-inv(M)*K,-inv(M)*C];
for i=1:N
    t1=t0+h;
    P=[0;0;0;0;0;0;0;0;0;inv(M)*[0;0;0;0;0;0;P0+PA*sin(w*t0);0;0]];
    k_1=f(t0,Y0,P);%函数定义参见f.m
    P=[0;0;0;0;0;0;0;0;0;inv(M)*[0;0;0;0;0;0;P0+PA*sin(w*(t0+h/2));0;0]];
    k_2=f(t0+h/2,Y0+h*k_1/2,P);
    k_3=f(t0+h/2,Y0+h*k_2/2,P);
    P=[0;0;0;0;0;0;0;0;0;inv(M)*[0;0;0;0;0;0;P0+PA*sin(w*(t0+h));0;0]];
    k_4=f(t0+h,Y0+h*k_3,P);
    Y1=Y0+(h/6)*(k_1+2*k_2+2*k_3+k_4);
    YY1_1(j)=Y1(1);YY1_2(j)=Y1(10);
    YY2_1(j)=Y1(2);YY2_2(j)=Y1(11);
    YY3_1(j)=Y1(3);YY3_2(j)=Y1(12);
    YY4_1(j)=Y1(4);
    YY4_2(j)=Y1(13);
    YY5_1(j)=Y1(5);YY5_2(j)=Y1(14);
    YY6_1(j)=Y1(6);YY6_2(j)=Y1(15);
    YY7_1(j)=Y1(7);YY7_2(j)=Y1(16);
    YY8_1(j)=Y1(8);YY8_2(j)=Y1(17);
    YY9_1(j)=Y1(9);YY9_2(j)=Y1(18);
    t0=t1;
    Y0=Y1;
    j=j+1;
end

参照这个，就可以完成你需要求解的微分方程啦！