Need Help Solving a system of 2nd and 1st Order ODEs

Question

Sherif Magdy 2020-12-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/702547-need-help-solving-a-system-of-2nd-and-1st-order-odes

评论： Sherif Magdy 2020-12-28

采纳的回答： Alan Stevens

在 MATLAB Online 中打开

I have three differenital equations modeling a suspension system

Variables are Xu, Xs, X3 while Xo is an input that is generated by a white noise. Everything else is a constant

I have been trying to solve it for some time now but for some reason the solutions are always either zero or completely wrong.

Any help would be greatly appreciated as I have only just started using MATLAB.

Here is my code:

function dy = mysys(t,y,T,Xo)
format long
Xo = interp1(T,Xo,t);
k1 = 15000; k2 = 15000; kt = 200000; c1 = 1000; c2 = 1000; Ms=300; Mu=40;
dy(1) = y(4);
dy(2) = y(5);
dy(3) = y(6);
dy(4) = ((k1/Ms)*(y(2)-y(3)) + (c1/Ms)*(y(5)-y(4)));
dy(3) = (y(5) - (k2/c2)*(y(3)-y(1)) - (k2/c2)*(y(2)-y(3)))
dy(6) = ((kt/Mu)*(Xo-y(2)) - (c1/Mu)*(y(5)-y(4)) - (c2/Mu)*(y(5)-y(6)) - (k2/Mu)*(y(2)-y(3)));
dy = dy(:);

timerange = [0 10];
y0 = [0 0 0 0 0 0];
T = road_xo_time.Time;
Xo = road_xo_time.Data;
[t,y] = ode45(@(t,y) mysys(t,y,T,Xo), timerange, init_cond);
plot(t, y(:,1))

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2020-12-27

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/702547-need-help-solving-a-system-of-2nd-and-1st-order-odes#answer_585092

在 MATLAB Online 中打开

I think you should be passing only 5, not 6, parameters to your mysys function, which is better written as

function dy = mysys(t,y,T,Xo)
    Xo = interp1(T,Xo,t);
    k1 = 15000; k2 = 15000; kt = 200000; c1 = 1000; c2 = 1000; Ms=300; Mu=40;
    xs = y(1);
    xu = y(2);
    x3 = y(3);
    vs = y(4);
    vu = y(5);
    v3 = vu - (k1/c2)*(x3 - xs) - (k2/c2)*(xu - x3);
    
    dy = [vs;           % dxs/dt
          vu;           % dxu/dt
          v3;           % dx3/dt
          (k1/Ms)*(xu - x3) + (c1/Ms)*(vu - vs);       % d2xs/dt2
          (kt/Mu)*(X0 - xu) - (c1/Mu)*(vu - vs) - (c2/Mu)*(vu - v3) - (k2/Mu)*(xu - x3)];   % d2xu/dt2
      
end