MATLAB: How to use ode45 to solve a system of two differential equation

mass spring systemode45second order odestwo degrees of freedom

I am aware of how to solve a system with one set of ODEs but not two

m1x1''=k2(x2-x1) – k1x1,

m2x2''=-k2(x2-x1) -k3x2

Best Answer

Feeling lazy today, and watching the Stuttgart-Hamburg football match, so I let the Symbolic Math Toolbox do the heavy lifting:

syms k1 k2 k3 m1 m2 x1(t) x2(t) x10 x20
Eqs = [m1*diff(x1(t),2) == k2*(x2(t)-x1(t)) - k1*x1(t),  m2*diff(x2(t),2) == -k2*(x2(t)-x1(t)) -k3*x2(t)];
SymSys = odeToVectorField(Eqs);
Sys = matlabFunction(SymSys)
Sys = @(Y,k1,k2,m1,m2)[Y(2);-(k2.*(Y(1)-Y(3)))./m2;Y(4);-(-k2.*Y(1)+k1.*Y(3)+k2.*Y(3))./m1];

It transmuted your ‘’x1 and ‘x2’ to elements of the ‘Y’ vector, but no worries. To make ode45 happy, you would change this to:

Sys = @(t,Y,k1,k2,m1,m2)[Y(2);-(k2.*(Y(1)-Y(3)))./m2;Y(4);-(-k2.*Y(1)+k1.*Y(3)+k2.*Y(3))./m1];

and then in ode45, call it as:

[t,y] = ode45(@(t,y) Sys(t,Y,k1,k2,m1,m2), tspan, Y0);

having previously defined the constants ‘k1’,‘k2’,‘m1’,‘m2’ in your workspace.

I did not run this, so I am labeling it UNTESTED CODE, but it should work.

Related Solutions

MATLAB: Use ODE45 to solve a system of two coupled second order ODEs

I prefer to let the Symbolic Math Toolbox do these derivations:

syms c1 c2 k1 k2 m1 m2 t Ftfcn x1(t) x2(t) Y 
%x1''=(F(t)-(c1+c2)*x1'+c2*x2'-(k1+k2)*x1+k2*x2)/m1
Dx1 = diff(x1);
D2x1 = diff(x1,2);
Dx2 = diff(x2);
D2x2 = diff(x2,2);
Eq1 = D2x1 == (Ftfcn-(c1+c2)*Dx1+c2*Dx2-(k1+k2)*x1+k2*x2)/m1
% x2''=(c2*x1'-c2*x2'+k2*x1-k2*x2)/m2
Eq2 = D2x2 == (c2*Dx1-c2*Dx2+k2*x1-k2*x2)/m2
[VF,Subs] = odeToVectorField(Eq1, Eq2)
ftotal = matlabFunction(VF, 'Vars',{t,Y,Ftfcn,c1,c2,k1,k2,m1,m2})

This produces (lightly edited):

ftotal = @(t,Y,Ftfcn,c1,c2,k1,k2,m1,m2)[Y(2);-(c2.*Y(2)-c2.*Y(4)+k2.*Y(1)-k2.*Y(3))./m2;Y(4);(Ftfcn(t)-(c1+c2).*Y(4)-(k1+k2).*Y(3)+c2.*Y(2)+k2.*Y(1))./m1];

Substituting random values and a random function:

Ftfcn = @(t) sin(t);
c1 = rand;
c2 = rand;
k1 = rand;
k2 = rand;
m1 = rand;
m2 = rand;
ic = zeros(4,1);
tspan = [0,1];
[T,Y] = ode45(@(t,Y) ftotal(t,Y,Ftfcn,c1,c2,k1,k2,m1,m2), tspan, ic);
plot(T, Y)
grid

This runs.

Experiment to get the result you want.

MATLAB: Solving a system of second order differential equations (coupled mass-spring system)

syms X1(t) X2(t) K f m1 m2 Y t
EqX1 = diff(X1,2) == K*(X1-X2)/m1;
EqX2 = diff(X2,2) == f + K*(X2-X1)/m2;
[ODE,Vars] = odeToVectorField(EqX1, EqX2)
Sys = matlabFunction(ODE, 'Vars', {t, Y, K, f, m1, m2})
m1 = 1;                 % Mass 1
m2 = 2;                 % Mass 2
K = 100;                % Spring constant
f = 5;                  % Forcing term on mass 2
v2_i = 3;
ic = [0 v2_i 0 0]';
tspan = [0, 1];
[T,Y] = ode45(@(t,Y)Sys(t,Y,K,f,m1,m2), tspan, ic);
plot(T,Y)

Best Answer

Related Solutions

MATLAB: Use ODE45 to solve a system of two coupled second order ODEs

MATLAB: Solving a system of second order differential equations (coupled mass-spring system)

Related Question