weixin_45266058 2021-12-19 17:03 采纳率: 100%
浏览 228
已结题

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

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

img

源程序如下:
%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

img


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

  • 写回答

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

    结果

    img

    本回答被题主选为最佳回答 , 对您是否有帮助呢?
    评论
查看更多回答(1条)

报告相同问题?

问题事件

  • 系统已结题 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广播数据遇到阻塞