How to use ode45 to solve a system of two differential equation?

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

Answer 1

[Star Strider]

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.

[Conner Nixon]

function Q3

syms k1 k2 k3 m1 m2 x1(t) x2(t) x10 x20

k1=25;

k2=25;

k3=25;

m1=1;

m2=1;

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);

tspan=[0 15];

Y0=[0 2];

Sys = matlabFunction(SymSys);

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];

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

end

This is my input, but I am unable to get it to run, could you tell me why this is so?

[Star Strider]

There are two problems, one mine (a typo in the ode45 call, the ‘@(t,y)’ should be ‘@(t,Y)’), the second that there need to be 4 initial conditions, since the ‘Sys’ function returns a (4x1) vector.

This works:

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];

k1=25;
k2=25;
k3=25;
m1=1;
m2=1;
tspan=[0 15];
Y0=[0 2 0 2];
[t,y] = ode45(@(t,Y) Sys(t,Y,k1,k2,m1,m2), tspan, Y0);

figure(1)
plot(t, y)
grid

You may want to change the initial conditions. I chose the second two arbitrarily.

[Conner Nixon]

Appologies, as I am not confident with Matlab as we have only been given limited teaching with it, but is it possible to isolate the plots of only Y(1) and Y(3)? Thank you for being so helpful!

[Star Strider]

My pleasure!

No apologies necessary. We’re all always learning.

To plot only ‘Y(1)’ and ‘Y(3)’ on one plot:

figure(1)
plot(t, Y(:,1))
hold on
plot(t, Y(:,3))
hold off
grid

The solutions are returned as column vectors, so it’s necessary to address them as columns of ‘Y’.

You could plot them using only one plot call and not using the hold function, but plotting them as I did here allows for more coding flexibility. I added a grid call because I like the reference lines.