weixin_45266058 2021-12-19 17:03 采纳率: 100%

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

%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)=25000
y(4) - 25000
y(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)=(307
y(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

• 写回答

#### 2条回答默认 最新

• 技术专家团-Joel 2021-12-19 20:24
关注

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

``````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))

``````

结果

本回答被题主选为最佳回答 , 对您是否有帮助呢?
评论

• 系统已结题 12月28日
• 已采纳回答 12月20日
• 修改了问题 12月19日
• 修改了问题 12月19日
• 展开全部

#### 悬赏问题

• ¥15 Sharepoint JS开发 付费技术指导
• ¥15 输入程序运行仿真后，烟雾值不实时检测，变成固定值
• ¥20 数据排序，可选择排序方向
• ¥15 修改一下代码，考虑进程到达时间不同的情况
• ¥15 帮我看看这是个啥题，带解题过程和结果，条件如下FCF = 290471.33 g1 = 15% r = 8% g2 = 4% n = 5
• ¥15 edem模拟颗粒不显示或者生成失败
• ¥15 Python代码编写
• ¥15 php 将rtmp协议转hls协议，无法播放
• ¥20 python代码编写
• ¥20 使用MPI广播数据遇到阻塞