Solving coupled second order differential equations
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Abraham Hartley
el 6 de Sept. de 2022
Comentada: Alan Stevens
el 6 de Sept. de 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 comentario
Torsten
el 6 de Sept. de 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.
Respuesta aceptada
Alan Stevens
el 6 de Sept. de 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 comentario
Alan Stevens
el 6 de Sept. de 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.
Más respuestas (0)
Ver también
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!