Plotting the derivative of infected population SI model

Hello;

I'm trying to draw the following model

I tried the code ode45 but it didn't work.

Is there a specific way to link to the same image? Thank you.

[Sam Chak]

@Amal Matrafi, could you please provide the ode45 code so that we can check? We are not experts in the SI model. Is 'SI' an abbreviation for something?

By the way, could you solve this ODE and plot the derivative of ?

[Amal Matrafi]

clc;
clear all
close all;
N=1000;
tend=300;
I0=10;
tspan = [0,tend];S0 = N - I0;y0 = [S0; I0];
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);

beta=0.1;
[t,y] = ode45(@(t,y) SIRfunc(t,y,beta,N), tspan, y0,opts);
plot(t,y(:,2),'b','LineWidth',1);
hold on

beta=0.2;
[t,y] = ode45(@(t,y) SIRfunc(t,y,beta,N), tspan, y0,opts);
plot(t,y(:,2),'r','LineWidth',1);
hold on

beta=0.25;
[t,y] = ode45(@(t,y) SIRfunc(t,y,beta,N), tspan, y0,opts);
plot(t,y(:,2),'LineWidth',1);
hold off
legend('\beta_{1}=0.1','\beta_{1}=0.2','\beta_{1}=0.25')

xlim([0 140])
ylim([0 1000])

function dydt = SIRfunc(~,y,beta,N)
dydt = [-beta/N*y(2)*y(1);
beta/N*y(2)*y(1)];
end

[Amal Matrafi]

My attempt

Answer 1

[Star Strider]

You need to plot the derivatives, not the solved equations.

One option is to use the gradient function:

dy2dt = gradient(y(:,2),t)

There are other, more direct (and probably more accurate) ways of calculating it from the original differential equation function. That requires a loop.

.

Answer 2

[Sam Chak]

Hi @Amal Matrafi

You can use deval() to obtain the first derivative. Alternatively, as suggested by @Star Strider, you can also use the gradient() approach to obtain the first derivative.

beta    = [0.1, 0.2, 0.25];
N       = 1000;
tend    = 150;
I0      = 10;
tspan   = [0,tend];
S0      = N - I0;
y0      = [S0; I0];
opts    = odeset('RelTol', 1e-2, 'AbsTol', 1e-4);

for j = 1:numel(beta)
    sol     = ode45(@(t, y) SIRfunc(t, y, beta(j), N), tspan, y0, opts);
    t       = linspace(0, 150, 1501);
    [y, yp] = deval(sol, t);
    plot(t, yp(2,:)), hold on
end
grid on
hold off
xlabel('t'), ylabel('dI/dt')
legend('\beta_{1} = 0.1','\beta_{2} = 0.2', '\beta_{3} = 0.25')

%% SI Model
function dydt = SIRfunc(t, y, beta, N)
    dydt = [-beta/N*y(2)*y(1);
             beta/N*y(2)*y(1)];
end

[Sam Chak]

Before learning to use deval(), I utilized the right-hand side of the state equation by directly substituting the solution from ode45(). This is pure math stuff!

beta    = [0.1, 0.2, 0.25];
N       = 1000;
tend    = 150;
I0      = 10;
tspan   = linspace(0, tend, 10*tend+1);
S0      = N - I0;
y0      = [S0; I0];
opts    = odeset('RelTol', 1e-2, 'AbsTol', 1e-4);

for j = 1:numel(beta)
    [t, y]  = ode45(@(t, y) SIRfunc(t, y, beta(j), N), tspan, y0, opts);
    dIdt    = beta(j)/N*y(:,2).*y(:,1);
    plot(t, dIdt), hold on
end
grid on
hold off
xlabel('t'), ylabel('dI/dt')
legend('\beta_{1} = 0.1','\beta_{2} = 0.2', '\beta_{3} = 0.25')

%% SI Model
function dydt = SIRfunc(t, y, beta, N)
    dydt = [-beta/N*y(2)*y(1);
             beta/N*y(2)*y(1)];
end