关于#matlab#的问题：[t,y]=ode45(@func,tspan,y0)

现象是利用ode函数解含变量振动耦合方程后，计算结果不收敛
该问题是为了解决五自由度振动含变量耦合方程的

源程序如下：
%ode45求解
y0=[0 0 0 0 0 0 0 0 0 0];
tspan=[0,37];
[t,y]=ode45(@func,tspan,y0);
figure(1)
plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5))

%func函数
function Dy=func(t,y)
Dy=zeros(10,1);
p1=-(sin(10t)(pi*(cos(10t -(pit)/40)/4800000+cos(10t +(pit)/40)/4800000)+cos(10t -(pit)/40)/12000-cos(10t +(pit)/40)/12000))/(pi^2- 160000)-(cos(10t)(- pi*(sin(10t -(pit)/40)/4800000+sin(10t +(pit)/40)/4800000)-sin(10t -(pit)/40)/12000+sin(10t +(pit)/40)/12000))/(pi^2- 160000)sin(tpi/40);
p2=-(sin(10t)((pi*(cos((3pi)/20- 10t +(pit)/40)/4800000+cos(10t +(3pi)/20+(pit)/40)/4800000)+cos((3pi)/20- 10t +(pit)/40)/12000-cos(10t +(3pi)/20+(pit)/40)/12000)))/(pi^2- 160000)-(cos(10t)(sin((3pi)/20- 10t +(pit)/40)/12000+sin(10t +(3pi)/20+(pit)/40)/12000- pi*(sin(10t +(3pi)/20+(pit)/40)/4800000-sin((3pi)/20- 10t +(pit)/40)/4800000)))/(pi^2- 160000)+(sin((3pi)/20)cos(10t))/(6000(pi^2- 160000) )sin(tpi/40);
p3=-(sin(10*(t+3))(pi(cos(10*(t+3) -(pi*(t+3))/40)/4800000+cos(10*(t+3) +(pi*(t+3))/40)/4800000)+cos(10*(t+3) -(pi*(t+3))/40)/12000-cos(10*(t+3) +(pi*(t+3))/40)/12000))/(pi^2- 160000)-(cos(10*(t+3))(- pi(sin(10*(t+3) -(pi*(t+3))/40)/4800000+sin(10*(t+3) +(pi*(t+3))/40)/4800000)-sin(10*(t+3) -(pi*(t+3))/40)/12000+sin(10*(t+3) +(pi*(t+3))/40)/12000))/(pi^2- 160000)sin((t+3)pi/40);
p4=-(sin(10(t+3))((pi*(cos((3pi)/20- 10(t+3) +(pi*(t+3)/40)/4800000+cos(10*(t+3) +(3pi)/20+(pi(t+3))/40)/4800000)+cos((3pi)/20- 10(t+3) +(pi*(t+3))/40)/12000-cos(10*(t+3) +(3pi)/20+(pi(t+3))/40)/12000)))/(pi^2- 160000)-(cos(10*(t+3))(sin((3pi)/20- 10*(t+3) +(pi*(t+3))/40)/12000+sin(10*(t+3) +(3pi)/20+(pi(t+3))/40)/12000- pi*(sin(10*(t+3) +(3pi)/20+(pi(t+3))/40)/4800000-sin((3pi)/20- 10(t+3) +(pi*(t+3))/40)/4800000)))/(pi^2- 160000)+(sin((3pi)/20)cos(10(t+3)))/(6000(pi^2- 160000) ))sin((t+3)pi/40);
Dy(1)=y(6);
Dy(2)=y(7);
Dy(3)=y(8);
Dy(4)=y(9);
Dy(5)=y(10);
Dy(6)=25000y(4) - 25000y(1) - 25000y(5) - (5y(6))/4 + (5y(9))/4 - (5y(10))/4;
Dy(7)=2500p1 + 2500p2 - 310y(2) + 600y(4) - 900y(5) - (35y(7))/2 + 175y(9) - (105y(10))/4;
Dy(8)=1250p3 + 1250p4 - (425y(3))/2 + (425y(4))/2 - (1275y(5))/4 - (425y(8))/2 + (425y(9))/2 - (1275y(10))/4;
Dy(9)=(200y(1))/3 + 16y(2) + (17y(3))/3 - (53y(4))/6 + (307y(5))/6 + y(6)/300 + (7y(7))/15 + (17y(8))/3 - (1841y(9))/300 - (2341y(10))/300;
Dy(10)=(307y(4))/16000 - (9y(2))/100 - (51y(3))/1600 - y(1)/4 - (521y(5))/3200 - y(6)/80000 - (21y(7))/8000 - (51y(8))/1600 - (2341y(9))/80000 + (4141*y(10))/80000;

end

希望有大神可以帮忙看一下是什么问题

解的性态就是这样吧，我用刚性求解一样是这个样子

M = diag(1./[80, 800,800,30000,8000000]);
C = [100, 0,0,-100,100;
    0,1400,0,-14000, 21000;
    0,0,170000,-170000,255000;
    -100,-14000,-170000,184100,234100;
    100,21000,255000,234100,-414100;
    ];
K = [
2000000,0,0,-2000000,2000000;
0,2480000,0,-480000,720000;
0,0,170000,-170000,255000;
-2000000,-480000,-170000,265000,-1535000;
2000000,720000,255000,-1535000,1302500
];
p1=@(t)-(sin(10*t)*(pi*(cos(10*t -(pi*t)/40)/4800000+cos(10*t +(pi*t)/40)/4800000)+...
    cos(10*t -(pi*t)/40)/12000-cos(10*t +(pi*t)/40)/12000))/(pi^2- 160000)-...
    (cos(10*t)*(- pi*(sin(10*t -(pi*t)/40)/4800000+sin(10*t +(pi*t)/40)/4800000)-...
    sin(10*t -(pi*t)/40)/12000+sin(10*t +(pi*t)/40)/12000))/(pi^2- 160000)*sin(t*pi/40);
p2=@(t)-(sin(10*t)*((pi*(cos((3*pi)/20- 10*t +(pi*t)/40)/4800000+...
    cos(10*t +(3*pi)/20+(pi*t)/40)/4800000)+cos((3*pi)/20- 10*t +...
    (pi*t)/40)/12000-cos(10*t +(3*pi)/20+(pi*t)/40)/12000)))/(pi^2- 160000)-...
    (cos(10*t)*(sin((3*pi)/20- 10*t +(pi*t)/40)/12000+sin(10*t +(3*pi)/20+(pi*t)/40)/12000-...
    pi*(sin(10*t +(3*pi)/20+(pi*t)/40)/4800000-sin((3*pi)/20- 10*t +(pi*t)/40)/4800000)))/(pi^2- 160000)+...
    (sin((3*pi)/20)*cos(10*t))/(6000*(pi^2- 160000) )*sin(t*pi/40);
p3=@(t)-(sin(10*(t+3))*(pi*(cos(10*(t+3) -(pi*(t+3))/40)/4800000+...
    cos(10*(t+3) +(pi*(t+3))/40)/4800000)+cos(10*(t+3) -...
    (pi*(t+3))/40)/12000-cos(10*(t+3) +(pi*(t+3))/40)/12000))/(pi^2- 160000)-...
    (cos(10*(t+3))*(- pi*(sin(10*(t+3) -(pi*(t+3))/40)/4800000+sin(10*(t+3) +...
    (pi*(t+3))/40)/4800000)-sin(10*(t+3) -(pi*(t+3))/40)/12000+sin(10*(t+3) +...
    (pi*(t+3))/40)/12000))/(pi^2- 160000)*sin((t+3)*pi/40);
p4= @(t)-(sin(10*(t+3))*((pi*(cos((3*pi)/20- 10*(t+3) +(pi*(t+3)/40)/4800000+cos(10*(t+3) +...
    (3*pi)/20+(pi*(t+3))/40)/4800000)+cos((3*pi)/20- 10*(t+3) +...
    (pi*(t+3))/40)/12000-cos(10*(t+3) +(3*pi)/20+(pi*(t+3))/40)/12000)))/(pi^2- 160000)-...
    (cos(10*(t+3))*(sin((3*pi)/20- 10*(t+3) +(pi*(t+3))/40)/12000+sin(10*(t+3) +...
    (3*pi)/20+(pi*(t+3))/40)/12000- pi*(sin(10*(t+3) +(3*pi)/20+...
    (pi*(t+3))/40)/4800000-sin((3*pi)/20- 10*(t+3) +...
    (pi*(t+3))/40)/4800000)))/(pi^2- 160000)+...
    (sin((3*pi)/20)*cos(10*(t+3)))/(6000*(pi^2- 160000) ))*sin((t+3)*pi/40);
Q = @(t)[
   0;
   2000000*(p1(t)+p2(t));
   1000000*(p3(t)+p4(t));
   0;
   0;
];
odefun = @(t, y) [
y(6:10);
M*(-C*y(6:10)-K*y(1:5)+Q(t));
];
%ode45求解
y0=[0 0 0 0 0 0 0 0 0 0]';
tspan=[0,37];
[t,y]=ode23s(odefun,tspan,y0);
figure(1)
plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5))

结果

大概率是这个系统发散。
将系统微分方程写成一阶状态方程形式，使用eig求出系统矩阵A的特征值，应该是有正实部的特征值。