清空筛选条件
清空筛选条件

Using a loop to get a script to repeat

2 次查看（过去 30 天）
Sam Hayes
Sam Hayes 2015-6-22
编辑： Sam Hayes 2015-6-22
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 个评论
Guillaume
Guillaume 2015-6-22
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Tim
Tim 2015-6-22
编辑：Tim 2015-6-22
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

类别

在 Help Center 和 File Exchange 中查找有关 Least Squares 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by