Solving coupled second order differential equations
4 views (last 30 days)
Show older comments
Abraham Hartley
on 6 Sep 2022
Commented: Alan Stevens
on 6 Sep 2022
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 Comment
Torsten
on 6 Sep 2022
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.
Accepted Answer
Alan Stevens
on 6 Sep 2022
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 Comment
Alan Stevens
on 6 Sep 2022
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.
More Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!