10自由度二阶微分方程组的MATLAB求解问题
我正在分析一个振动问题，列出的是10自由度二阶微分方程组，形如：{X''}+{X'}+{X}={P}
其中=diag{J1,J2,J3,J4,J5,J6,J7,J8,m,m}
{X}=（φ1，φ2，φ3，φ4，φ5，φ6，φ7，φ8，X,Y）
为10*10阶的常数矩阵
为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
回复 #2 花如月 的帖子
天哪，十自由度、二阶微分方程组，兄弟，你比我还狠，这个太复杂了，解方程不好解啊！你要是用ode方法的话，就参考本板块的一个关于微分方程的精华贴，照着做就可以了，不过解这样的方程可是需要时间的啊！要耐心哦！ 将原问题贴出来看看.
可以先试试ode45之类的函数． 下面的程序是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=;
K=;
C=zeros(9,9);
%C=0.01*;
I=eye(9,9);
A=;
ydot=A*Y+P; 第二部分：求解
clear;
clc;
t0=0;
Y0=;
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=;
K=;
%C=zeros(9,9);
C=0.01*;
I=eye(9,9);
A=;
for i=1:N
t1=t0+h;
P=];
k_1=f(t0,Y0,P);%函数定义参见f.m
P=];
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=];
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
参照这个，就可以完成你需要求解的微分方程啦！
页:
[1]