Solving coupled second order differential equations

12 ビュー (過去 30 日間)
Abraham Hartley
Abraham Hartley 2022 年 9 月 6 日
コメント済み: 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 件のコメント
Torsten
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.

サインインしてコメントする。

採用された回答

Alan Stevens
Alan Stevens 2022 年 9 月 6 日
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 件のコメント
Alan Stevens
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.

サインインしてコメントする。

その他の回答 (0 件)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by