# Custom Integrator

We will use the same differential equations for the SIR model as previously introduced. The right hand side of the ordinary differential equations are evaluated with the following function eval_SIR_rhs.m:

```function xdot = eval_SIR_rhs(t, x, c)
% Inputs:
% t : 1x1 scalar time
% x : 3x1 state vector [S; I; R]
% c : structure containing 2 constant parameters
% Outputs:
% xdot : 3x1 state derivative vector [S'; I', R']

S = x(1); % susceptible [person]
I = x(2); % infected [person]
R = x(3); % recovered [person]

beta = c.beta;  % infectious contacts per time [person/day]
tau = c.tau;  % mean infectious period [day]

N = S + I + R;

Sdot = -beta*I*S/N;

Idot = beta*I*S/N - I/tau;

Rdot = I/tau;

xdot = [Sdot; Idot; Rdot];

end
```

The following function euler_integrate.m implements Euler's Method of integration and is designed to be a drop in replacement for ode45():

```function [ts, xs] = euler_integrate(f, ts, x0)
% f: function handle for the state derivative function
% ts: 1xn vector of time values
% x0: mx1 vector of initial states
% xs: nxm matrix where each column is the state at the ith time

xs = zeros(length(ts), length(x0));  % nxm matrix of zeros

xs(1, :) = x0;  % set first row to initial condition

for i=2:length(ts)
deltat = ts(i) - ts(i - 1);
% note the transpose ' operators to ensure the dimensions match,
% as xs(i - 1, :) is 1xm instead of mx1
xs(i, :) = xs(i - 1, :)' + deltat * f(ts(i - 1), xs(i - 1, :)');
end

end
```

The ordinary differential equations are integrated (simulated) with the following script simulate_SIR_model_euler.m:

```% constant parameters
c.beta = 1/7;  % infectious contacts per day [person/days]
c.tau = 14;  % days infectious [day]

% intial population values
S0 = 39.56E6;  % population of California [person]
I0 = 1000;  % infected [person]
R0 = 200;  % recovered [person]

x0 = [S0; I0; R0];

num_months = 8;
num_days_per_month = 30;

ts = 0:num_months*num_days_per_month;  % units of days

% integrate the equations with both the Euler method and ode45 (Runga-Kutta
% method)
f = @(t, x) eval_SIR_rhs(t, x, c);
[ts, xs] = euler_integrate(f, ts, x0);

[ts2, xs2] = ode45(f, ts, x0);

% plot the results, comparing the two integration results
hold on
title(sprintf('R0 = %1.2f', c.beta*c.tau))
plot(ts, xs(:, 1), 'b', ts2, xs2(:, 1), 'b--')  % S(t)
plot(ts, xs(:, 2), 'r', ts2, xs2(:, 2), 'r--')  % I(t)
plot(ts, xs(:, 3), 'g', ts2, xs2(:, 3), 'g--')  % R(t)
xlabel('Time [days]')
ylabel('Number of people')
legend('S (Euler)', 'S (ode45)', ...
'I (Euler)', 'I (ode45)', ...
'R (Euler)', 'R (ode45)')
hold off
```

And produces the following figure: The Euler Method deviates from the ode45() solution due to the poorer error properties of the method.

The following was not covered in lab, but you may find it interesting.

A simple Runga-Kutta Method of order 4 can also be implemented and compared to the Runga-Kutta method used by ode45(). The file runga_kutta_integrate.m is shown below:

```function [t, x] = runga_kutta_integrate(f, t, x0)
% RUNGA_KUTTA_INTEGRATE - Integrates a set of first order ordinary
% differential equations using the 4th order Runga-Kutta method.
%
% Syntax: [t, x] = runga_kutta_integrate(f, t, x0)
%
% Inputs:
%   f - An anonymous function that evaluates the right hand side of the
%       first order ordinary differential equations, i.e. dx/dt. The
%       function must follow the syntax dxdt = f(t, x) and return a size
%       mx1 vector where t is size 1x1 and x is size mx1.
%   t - Vector of monotonically increasing time values. [size 1xn]
%   x0 - Vector a initial conditions for the states. [size mx1]
% Outputs:
%   t - Vector of monotonically increasing time values (same as the
%       input t). [size 1xn]
%   x - Matrix of state values as a function of time. [size nxm]

% Create an "empty" matrix to hold the results n x m.
x = nan * ones(length(t), length(x0));

% Set the initial conditions to the first element.
x(1, :) = x0;

% Use a for loop to sequentially calculate each new x.
for i = 2:length(t)
deltat = t(i) - t(i-1);

k1 = f(t(i - 1), x(i - 1, :)');
k2 = f(t(i - 1) + deltat/2, x(i - 1, :)' + deltat.*k1./2);
k3 = f(t(i - 1) + deltat/2, x(i - 1, :)' + deltat.*k2./2);
k4 = f(t(i), x(i - 1, :)' + deltat.*k3);

x(i, :) = x(i - 1, :)' + 1/6*deltat*(k1 + 2*k2 + 2*k3 + k4);

end

end
```

The ordinary differential equations are integrated (simulated) with the following script simulate_SIR_model_ode45_euler.m:

```% constant parameters
c.beta = 1/7;  % infectious contacts per day [person/days]
c.tau = 14;  % days infectious [day]

% initial population values
S0 = 39.56E6;  % population of California [person]
I0 = 1000;  % infected [person]
R0 = 200;  % recovered [person]

x0 = [S0; I0; R0];

num_months = 8;
num_days_per_month = 30;

ts = 0:num_months*num_days_per_month;  % units of days

% integrate the equations
f = @(t, x) eval_SIR_rhs(t, x, c);
[ts, xs] = ode45(f, ts, x0);

[ts2, xs2] = euler_integrate(f, ts, x0);

[ts3, xs3] = runga_kutta_integrate(f, ts, x0);

% plot the results
hold on
title(sprintf('R0 = %1.2f', c.beta*c.tau))
plot(ts, xs(:, 1), 'b', ts2, xs2(:, 1), 'b--', ts3, xs3(:, 1), 'b.')  % S(t)
plot(ts, xs(:, 2), 'r', ts2, xs2(:, 2), 'r--', ts3, xs3(:, 2), 'r.')  % I(t)
plot(ts, xs(:, 3), 'g', ts2, xs2(:, 3), 'g--', ts3, xs3(:, 3), 'g.')  % R(t)
xlabel('Time [days]')
ylabel('Number of people')
legend('Susceptible (ode45)', 'Susceptible (Euler)', 'Susceptible (RK)', ...
'Infected (ode45)', 'Infected (Euler)', 'Infected (RK)', ...
'Recovered (ode45)', 'Recovered (Euler)', 'Recovered (RK)')
hold off
```

And produces the following figure: The custom Runga-Kutta method is essentially identical to the ode45() result for this system.