[Math] Solving a system of ODE (3 equations) using runge-kutta method order 4 – Matlab

MATLAB

$$x' = \lambda – \rho x – \beta xz \\
y' = \beta xz – \delta y \\
z' = py – cz$$

$x_0=43100, y_0 = 0, z_0 = 0.0033$

$\lambda = 388, \rho = 0.009, \beta =3.61e-8, \delta = 0.18, p = 50000, c = 23. $

my code is as follows:

f = @(t,x) [388-0.009*x(1)-0.0000000361*x(1)*x(3); 0.0000000361*x(1)*x(3) - 0.18*x(2); 50000*x(3) - 23*x(3)];
tspan = 0:7:84;
[t,xa] = ode45(f,tspan,[43100 0 0.0033]);

Issue:
The code keeps running without any result or termination. So I have tried using ode23s and ode15s but I get a number of warning messages like

In ode23s (line 379) Warning: Matrix is close to singular or badly
scaled. Results may be inaccurate. RCOND = 7.484584e-246. Warning:
Failure at t=7.014056e+00. Unable to meet integration tolerances
without reducing the step size below the smallest value allowed
(2.491893e-14) at time t.

I am wondering if my code is incorrect? Please help!

Best Answer

After trying to re-scale your system I've just found a little typo in your code in a last component of f(t,x). Now everything is working, but if you want you can use my version of your code, it looks different, but after all all the changes just minor.

lambda = 388;
rho = 0.009;
beta=3.61e-8;
delta=0.18;
p=50000;
c=23;

f = @(t,x) [lambda-rho*x(1)-beta*x(1)*x(3); beta*x(1)*x(3) - delta*x(2); p*x(2) - c*x(3)];
tspan = [0,1];
options = odeset('RelTol',1e-1,'AbsTol',[1e-4 1e-4 1e-4]);
[t,xa] = ode45(f,tspan,[43100 0 0.0033],options);
plot(t,xa(:,1),'-'); figure
plot(t,xa(:,2),'-.');figure
plot(t,xa(:,3),'.');

Related Solutions

[Math] Using Runge Kutta 2 on a System of Equations

Note that I have not checked you have implemented the RK2 algorithm correctly. I have just modified the code that you have provided.

You don't need to change very much at all.

You should preallocate (which is good practice) y = zeros(2,nstep+1);, where column n of y is the solution at time n-1.

You need to set the initial condition:

y(:,1) = [value1;value2];

Your expression for k1 needs to change so that the vector of current y values is passed to f:

k1 = h * fxn(t(n), y(:,n));

Your expression for k2 needs to change similarly:

k2 = h * fxn(t(n) + h/2, y(:,n) + k1/2);

And you need to update when you store the solution

y(:,n+1) = y(:,n) + k2;

This will store the solution for y1 in the first row of y and y2 in the second row.

As you store the solution for each timestep as a column vector, you need to ensure that your function fxn returns a column vector.

[Math] Runge Kutta Method Matlab code

You are exactly correct in your assumption. You should make this into a system of first-order ODEs. The straight forward way to transform this equation into a first order ODE is to create the vector $X(t) \in \mathbb{R}^{n}$ where each element is defined as $X(i) := x^{(i)}(t)$.

Then you have a first-order equation defined as

$$X' = M X$$

where $M \in \mathbb{R}^{nxn}$. You can fully describe $M$ noting that for every $i$, by construction:

$$(X(i-1))' = (x^{(i-1)}(t))' = x^{(i)}(t) = X'(i)$$

This defines $n-1$ equations, with the last equation in $M$ being your original ODE re-written in terms of $X$. Note also that you have $X(t_0)$ as an input. This is all the information you need to solve your ODE using RK4.

Best Answer

Related Solutions

[Math] Using Runge Kutta 2 on a System of Equations

[Math] Runge Kutta Method Matlab code

Related Question