How to solve SIR model with using DTM (Differential Transform Method)

11 views (last 30 days)
I am trying to use dtm for solving SIR model, Although my code is run but I think the DTM part is wrong. I need help for DTM part
here is my code
function sir_model_dtm()
% Parameters
alpha = 0.3; % Infection rate
beta = 0.1; % Recovery rate
N = 1000; % Total population
I0 = 1; % Initial number of infected individuals
R0 = 0; % Initial number of recovered individuals
S0 = N - I0 - R0; % Initial number of susceptible individuals
% Time parameters
T = 100; % Total simulation time
dt = 1; % Time step
% Initialize arrays to store results
S = zeros(1, T);
I = zeros(1, T);
R = zeros(1, T);
% Initial conditions
S(1) = S0;
I(1) = I0;
R(1) = R0;
% Simulate the SIR model using DTM
for t = 2:T
% Compute new values
S(t) = S(t-1) - alpha*S(t-1)*I(t-1)/N * dt;
I(t) = I(t-1) + (alpha*S(t-1)*I(t-1)/N - beta*I(t-1)) * dt;
R(t) = R(t-1) + beta*I(t-1) * dt;
end
% Plot the results
t = 0:dt:T-dt;
plot(t, S, 'b', t, I, 'r', t, R, 'g', 'LineWidth', 2);
legend('Susceptible', 'Infectious', 'Recovered');
xlabel('Time');
ylabel('Number of individuals');
title('SIR Model');
end

Answers (1)

Athanasios Paraskevopoulos
function sir_model_dtm()
% Parameters
alpha = 0.3; % Infection rate
beta = 0.1; % Recovery rate
N = 1000; % Total population
I0 = 1; % Initial number of infected individuals
R0 = 0; % Initial number of recovered individuals
S0 = N - I0 - R0; % Initial number of susceptible individuals
% Time parameters
T = 100; % Total simulation time
dt = 1; % Time step
% Initialize arrays to store results
S = zeros(1, T);
I = zeros(1, T);
R = zeros(1, T);
% Initial conditions
S(1) = S0;
I(1) = I0;
R(1) = R0;
% Simulate the SIR model using DTM
for t = 2:T
% Current values
S_curr = S(t-1);
I_curr = I(t-1);
R_curr = R(t-1);
% Compute intermediate values using DTM (predictor step)
S_inter = S_curr - alpha * S_curr * I_curr / N * dt;
I_inter = I_curr + (alpha * S_curr * I_curr / N - beta * I_curr) * dt;
R_inter = R_curr + beta * I_curr * dt;
% Compute new values (corrector step)
S(t) = S_curr - 0.5 * dt * (alpha * S_curr * I_curr / N + alpha * S_inter * I_inter / N);
I(t) = I_curr + 0.5 * dt * ((alpha * S_curr * I_curr / N - beta * I_curr) + (alpha * S_inter * I_inter / N - beta * I_inter));
R(t) = R_curr + 0.5 * dt * (beta * I_curr + beta * I_inter);
end
% Plot the results
t = 0:dt:T-dt;
plot(t, S, 'b', t, I, 'r', t, R, 'g', 'LineWidth', 2);
legend('Susceptible', 'Infectious', 'Recovered');
xlabel('Time');
ylabel('Number of individuals');
title('SIR Model');
end

Categories

Find more on Optimization Toolbox in Help Center and File Exchange

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!