麻烦大家帮忙看看这个振动微分方程，不知道如何去解......

冬夜

     大家好，我的这个振动微分微分方程里面由于质量是耦合的，在用ode求解时，由于赋初值时只有y0，再代入方程迭代时，就算不了，不知道ode是不不能解这类方程，还请大师指导！谢谢大家，以下是方程及程序代码，期待得到大家的指导，真心求教......
     函数程序：
function dy=xuanjiafun(t,y)
m=1380;I=2440;J=380;e=0;i=0;l1=1.25;l2=1.51;
b1=0.74;b2=0.74;k11=17;k21=17;k31=17;k41=17;c1=1.5;c2=1.5;c3=1.5;c4=1.5;A=0.05;w=8;t1=0.13;
dy(1)=y(2);
dy(3)=y(4);
dy(5)=y(6);
dy(2)=(-(k11*(y(1)-l1*y(3)-b1*y(5))+k21*(y(1)-l1*y(3)+b2*y(5))+k31*(y(1)+l2*y(3)-b1*y(5))+k41*(y(1)+l2*y(3)+b2*y(5))+c1*(dy(1)-l1*dy(3)-b1*dy(5))+c2*(dy(1)-l1*dy(3)+b2*dy(5))+c3*(dy(1)+l2*dy(3)-b1*dy(5))+c4*(dy(1)+l2*dy(3)+b2*dy(5)))+2*A*sin(w*t)+2*A*sin(w*(t+t1)))/m-e*dy(4)-i*dy(6);
dy(4)=((-((b1^2+l1^2)^1/2+(b2^2+l1^2)^1/2)*A*sin(w*t)+((b2^2+l2^2)^1/2+(b1^2+l2^2)^1/2)*A*sin(w*(t+t1))-(l2*(k31*(y(1)+l2*y(3)-b1*y(5))+k41*(y(1)+l2*y(3)+b2*y(5)))-l1*(k11*(y(1)-l1*y(3)-b1*y(5))+k21*(y(1)-l1*y(3)+b2*y(5)))-l1*(c1*(dy(1)-l1*dy(3)-b1*dy(5))+c2*(dy(1)-l1*dy(3)+b2*dy(5)))+l2*(c3*(dy(1)+l2*dy(3)-b1*dy(5))+c4*(dy(1)+l2*dy(3)+b2*dy(5))))-I*dy(4))/(e*m)-dy(2)-i*dy(6))/e;
dy(6)=(((b2^2+l1^2)^1/2-(b1^2+l1^2)^1/2)*A*sin(w*t)+((b2^2+l2^2)^1/2-(b1^2+l2^2)^1/2)*A*sin(w*(t+t1))-(b2*(k21*(y(1)-l1*y(3)+b2*y(5))+k41*(y(1)+l2*y(3)+b2*y(5)))-b1*((k11*(y(1)-l1*y(3)-b1*y(5))+k31*(y(1)+l2*y(3)-b1*y(5)))-b1*(c1*(dy(1)-l1*dy(3)-b1*dy(5))+c3*(dy(1)+l2*dy(3)-b1*dy(5)))+b2*(c2*(dy(1)-l1*dy(3)+b2*dy(5))+c4*(dy(1)+l2*dy(3)+b2*dy(5))))-J*dy(6))/(i*m)-dy(2)+e*dy(4))/i;
t(end,1)
dy=[dy(1);dy(2);dy(3);dy(4);dy(5);dy(6)];
主程序：
clc
clear
hold on;
box on;
Y1=[];
w=8;k=0;w_vc =w/(2*pi);period=1/w_vc;
step=period/100;
tspan=0:step:500*period;
y0=[0;0;0;0;0;0];
options=odeset('RelTol',1.0e-9); % 10-9 for ode15s; 10-5 for ode45  
[t,y]=ode45('xuanjiafun',tspan,y0,options);   %计算积分
    %[t,y]=ode45('bevelchaos_fun',tspan,y0,options,b);   %计算积分
Y100=[y(25000:1000:end,1)];
Y200=[y(25000:1000:end,2)];
Y300=[y(25000:1000:end,3)];
Y400=[y(25000:1000:end,4)];
Y500=[y(25000:1000:end,5)];
Y600=[y(25000:1000:end,6)];

save Y100.dat Y100 -ascii
save Y200.dat Y200 -ascii
save Y300.dat Y300 -ascii
save Y400.dat Y400 -ascii
save Y500.dat Y500 -ascii
save Y600.dat Y600 -ascii
y=[t,y]
save datay.dat y -ascii

飞天小男警

同求同求！质量耦合还能不能用龙格库塔法？？求大神

dw04116

希望你用公式说明，或者文字。。看程序很累

冬夜

dw04116 发表于 2013-3-13 11:39

                               
登录/注册后可看大图

希望你用公式说明，或者文字。。看程序很累

我将没用的删掉，方程如下：
function dy=fun(t,y)
dy(1)=y(2);
dy(3)=y(4);
dy(5)=y(6);
dy(2)=。。。。。。-e*dy(4)-i*dy(6);
dy(4)=。。。。。。。。
dy(6)=。。。。。。。。
t(end,1)
dy=[dy(1);dy(2);dy(3);dy(4);dy(5);dy(6)]
在用龙格库塔计算时，赋初值是y0=[0;0;0;0;0;0]，在算dy(2)的时候出现了dy(4)和dy(6)由于没有初值算不下去，还是我在龙格库塔解法的解法上理解有问题？

冬夜

问题解决啦！小小的问题......还是要认真学习，哈哈......

godlinux

冬夜 发表于 2013-3-14 11:03

                               
登录/注册后可看大图

问题解决啦！小小的问题......还是要认真学习，哈哈......

怎么解决的？求指导！谢谢