Using a loop to get a script to repeat
2 views (last 30 days)
Show older comments
Hi, I'm modelling disease spread using the SIR model and I've managed to code an m.file that works, but now I need to be able to get the model to loop, say 100 times, changing the value of the parameter beta by 0.01 each time to then get a range of values for R_0 which I can then put in a vector and plot to hopefully see a bifurcation.
Is there a way I can do this using loops rather than brute force going through each value of beta and then writing down the corresponding value of R_0?
This is my code:
function SIA
param.beta = 0.3; % force of infection
param.mu = 0.5; % birth rate
param.nu = 0; % death rate
param.delta = 0.8; % virulence of AIDS
param.gamma = 0.5; % proportion who develop AIDS
param.alpha= 0; % virulence of HIV
% Initial conditions are the values of all variables at time zero.
initial.S = 10^6; % set the initial value of 'S'
initial.I = 10^5; % set the initial value of 'I'
initial.A = 10^3; % set the initial value of 'A'
inits = [initial.S; initial.I; initial.A]; % Initial conditions
end_time = 50; % set how long to run for
function deriv = ode_system (t, x, param);
% Function to calculate derivatives of the SIR model
% Define N and calculate R_0
N = initial.S + initial.A + initial.I; % Population size
R_0 = param.beta / (param.gamma+param.nu+param.alpha)
% Define ODEs
S = x(1); % Number of susceptibles
I = x(2); % Number of HIV infected
A = x(3); % Number of AIDS infected
dS = +param.mu*N-param.beta * S * I/N-param.nu*S;
dI = +param.beta * S * I/N - (param.nu+param.gamma+param.alpha) * I;
dA = +param.gamma * I-(param.nu+param.delta)*A;
deriv = [dS; dI; dA];
end
% integrate the ODE system
[t, y] = ode45(@(t, x) ode_system(t, x, param),[0 end_time],inits,[]);
1 Comment
Guillaume
on 22 Jun 2015
Please, remove all those blank lines that you probably spent a while putting in between each line of code and then click on the {} Code button to format them as code. Much faster and more readable this way.
Answers (1)
Tim
on 22 Jun 2015
Edited: Tim
on 22 Jun 2015
Change "function SIA" to accept an input and return the integration. then call it from another script like so:
b=1:.01:whatever;
for i = 1:100
[t,y]=SIA(b(i));
result(i,1)= t;
result(i,2)= y;
end
and SIA will now be:
function [t,y] = SIA(b)
param.beta = b;
blah
blah
blah
0 Comments
See Also
Categories
Find more on Least Squares in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!