How do you use ODE to solve four linear differential equations?

I have two coupled systems composed of four linear differential equations. How do you use a stiff ode to solve these four equations simultaneously?

m is a constant value

x1' = m * [y1 + x1 - x2];

y1' = -1/m * [x1 + y1];

x2' = m * [y2 + x2 - x1];

y2' = -1/m * [x2 + y2];

I solved the differential equation for one system, the code is as follows:

f = @(t,y) [m * (y(2) + y(1) - z); (-1/m) * ( y(1) + y(2)) ];

g = @(y) f(0,y);

[T Y] = ode15s(f, [0 80], [0, 0]);

In this case, z was a constant. However, the solution of the first system is now dependent on the second (x2).

Answer 1

[Star Strider]

I would take advantage of the Symbolic Math Toolbox to create an anonymous function you can then use directly with any of the numerical ODE solvers.

syms x1(t) x2(t) y1(t) y2(t) t Y m 
Eq1 = diff(x1) == m*(y1 + x1 + x2);
Eq2 = diff(y1) == -1/m*(x1 + y1);
Eq3 = diff(x2) == m*(y2 + x2 + x1);
Eq4 = diff(y2) == -1/m*(x2 + y2);
[DEVF, Subs] = odeToVectorField(Eq1, Eq2, Eq3, Eq4);
ODEfcn = matlabFunction(DEVF, 'Vars',{t,Y,m})

producing (after some slight editing):

ODEfcn = @(t,Y,m)[-(Y(1)+Y(2))./m;m.*(Y(1)+Y(2)+Y(3));m.*(Y(2)+Y(3)+Y(4));-(Y(3)+Y(4))./m];

You would then call this in the ODE solvers as:

[t,y] = ode15s(@(t,Y) ODEfcn(t,Y,M), tspan, ic);

where ‘tspan’ is the vector of times or time limits, and ‘ic’ is a (4x1) vector of initial conditions.

Note that the ‘Subs’ variable describes the substitutions made to create the rows of ‘DEVF’. I prefer using the numeric ODE solvers.

You can also use the dsolve (link) function if you only want an analytic solution. For best results, you have to give it the initial values for all the derivatives.

syms x1(t) x2(t) y1(t) y2(t) t Y m X10 X20 Y10 Y20
Eq1 = diff(x1) == m*(y1 + x1 + x2);
Eq2 = diff(y1) == -1/m*(x1 + y1);
Eq3 = diff(x2) == m*(y2 + x2 + x1);
Eq4 = diff(y2) == -1/m*(x2 + y2);
[x1,x2,y1,y2] = dsolve(Eq1, Eq2, Eq3, Eq4, x1(0) == X10,  x2(0) == X20, y1(0) == Y10, y2(0) == Y20)

[EL890]

That worked perfectly. Thanks!

[Star Strider]

As always, my pleasure!

[Mustafain Ali]

Hello,

Thank you for explaining the method to analytically solve an ODE system through the symbolic toolbox. However, I am trying to do the same as you did but I am getting an error:

syms y1(t) y2(t) t y1 y2 y10 y20

eqn1 = diff(y1) == -y1 - y1*y2^2 + 294*y2;

eqn2 = diff(y2) == (y1 - y1*y2)/98 - 3*y2;

[y1, y2] = dsolve(eqn1,eqn2, y1(0) == y10, y2(0) == y20)

It gives the error:

Array indices must be positive integers or logical values.

Error in indexing (line 968)

R_tilde = builtin('subsref',L_tilde,Idx);

Error in stiffODE (line 18)

[y1, y2] = dsolve(eqn1,eqn2, y1(0) == y10, y2(0) == y20)

Am I missing something?

[Star Strider]

@Mustafain Ali — A few things, actually.

The first is that the differential equations are nonlinear (the product makes them so), so it is highly unlikely they have an analytic solution (most nonlinear differential equations do not).

Beyond that, you first declared ‘y1(t)’ and ‘y2(t)’ as fuctions of time (appropriate) and then defined them again as scalars ‘y1’ and ‘y2’ (in error), and that threw the indexing error.

Correcting that (and eliminating the two initial conditions that are no longer relevant) you can then integrate the system numerically using the steps and functions I use here —

syms y1(t) y2(t) t Y
eqn1 = diff(y1) == -y1 - y1*y2^2 + 294*y2;
eqn2 = diff(y2) == (y1 - y1*y2)/98 - 3*y2;
[VF,Sbs] = odeToVectorField(eqn1, eqn2)

VF = 

Sbs = 

odefcn = matlabFunction(VF, 'Vars',{t,Y})

odefcn = function_handle with value:
    @(t,Y)[Y(1).*-3.0+Y(2)./9.8e+1-(Y(1).*Y(2))./9.8e+1;-Y(1).^2.*Y(2)+Y(1).*2.94e+2-Y(2)]

[t,y] = ode45(odefcn, [0 20], [1 1]);

figure
yyaxis left
plot(t, y(:,2))
ylabel('y_2(t)')
yyaxis right
plot(t, y(:,1))
ylabel('y_1(t)')
grid
xlabel('t')
legend(string(Sbs))

Choose appropriate initial conditions, then do the numeric integration.

.