Solving coupled second order differential equations

Question

Abraham Hartley 2022-9-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1795765-solving-coupled-second-order-differential-equations

评论： Alan Stevens 2022-9-6

Hi,

I constantly get an error that I do not know how to resolve when I run this code. Any help would be greatly appreciated. The Force function applies a force of 50N @ t = 0 and equals 0 at all other time steps.

syms c k m1 m2 x1(t) x2(t) t F Y;

% 1st and 2nd derivative

dx1 = diff(x1);

d2x1 = diff(x1,2);

dx2 = diff(x2);

d2x2 = diff(x2,2);

% Defining equations

Eq1 = d2x1 == -(c/m1)*(dx1-dx2) - (k/m1)*(x1(t)-x2(t)) + F/m1;

Eq2 = d2x2 == -(c/m2)*(dx1-dx2) + (k/m2)*(x1(t)-x2(t));

[VF,subs] = odeToVectorField(Eq1, Eq2);

ftotal = matlabFunction(VF,'Vars',{t,Y,F,c,k,m1,m2});

m1 = 10;

m2 = 20;

c = 30;

k = 5;

F = @(t) Force(t);

tspan = [0 1];

ic = [8 0 0 0];

[t,Y] = ode45(@(t,Y) ftotal(t,Y,F,c,k,m1,m2), tspan, ic);

figure

plot(t, Y)

grid

legend(string(subs))

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2022-9-6

The value of the "Force" function at t=0 must somehow be part of the initial conditions ic.

As part of the differential equations, it won't influence the result.

请先登录，再进行评论。

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

Answer 1

Alan Stevens 2022-9-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1795765-solving-coupled-second-order-differential-equations#answer_1042885

在 MATLAB Online 中打开

Something like this?

m1 = 10;

m2 = 20;

c = 30;

k = 5;

tspan = [0 1];

ic = [8 0 0 0];

[t,Y] = ode45(@(t,Y) ftotal(t,Y,c,k,m1,m2), tspan, ic);

figure

plot(t, Y)

xlabel('t'), ylabel('y and dydt')

grid

legend('y1','dy1/dt','y2','dy2dt')

function dYdt = ftotal(t,Y,c,k,m1,m2)

y1 = Y(1); dy1dt = Y(2);

y2 = Y(3); dy2dt = Y(4);

dYdt = [dy1dt;

-(c/m1)*(dy1dt - dy2dt) - (k/m1)*(y1 - y2) + 50*(t==0)/m1;

dy2dt;

-(c/m2)*(dy1dt - dy2dt) + (k/m2)*(y1 - y2)];

end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Alan Stevens 2022-9-6

Instead of 50*(t==0)/m1; you could try 50*(t<0.01)/m1; but you still won't see an effect. However, if you try 5000*(t<0.01)/m1; you start to see an effect. If you try decreasing the 5000 to smaller values, towards 50, you see the effect getting smaller and smaller.

请先登录，再进行评论。