Help with coding a model
1 回表示 (過去 30 日間)
古いコメントを表示
Greetings,
I am currently trying to simulate an SI (SIR) model but everytime I get a new error and I do not know how to solve it from now on, could any one please help me?
The model is simple as the following:

With beta = 0.007
gamma = 0.01
L = 80
And 

Sigma = 0.115129254649702
I tried using Euler method and ode45 but with no luck, I am a bit confused on how to discretize just over time and inside the integral I should descritize over space
Here is my current code
syms x
L=80;
y=[0 L];
j = [0,L];
days = 10;
tspan = [1 10];
gamma=10^-2;
R_0 = 0.7;
beta=R_0*gamma;
epsilon = 0.0001;
sigma = -log(epsilon)/L;
Z0 = 0.2;
S0 = 1 - Z0;
[t,Z] = ode45(@(t,Z) odefun(Z,x,y,L,sigma,beta,gamma,days),tspan,Z0);
grid on;
hold on
function dZdt = odefun(Z,x,y,L,sigma,beta,gamma,days)
dZdt = zeros(L,days);
dZdt = beta*(1-Z)*int(exp(-sigma*abs(x-y))*Z(y),y,0,L) - gamma*Z;
end
Thank you in advance!
Have a nice day
0 件のコメント
採用された回答
Davide Masiello
2022 年 10 月 9 日
編集済み: Davide Masiello
2022 年 10 月 9 日
I think this could be the solution to your problem
% Parameters
L = 80;
gamma = 1e-2;
R_0 = 0.7;
beta = R_0*gamma;
epsilon = 0.0001;
sigma = -log(epsilon)/L;
% Time discretization
days = 10;
tspan = [1 days];
% Space discretization
y = linspace(0,L,100);
% Intial conditions
z0 = 0.2*ones(size(y));
[t,Z] = ode45(@(t,z) odefun(t,z,y,sigma,beta,gamma),tspan,z0);
plot(y,Z) % plot of Z vs space at various times
xlabel('y')
ylabel('Z')
legend([repmat('t = ',length(t),1),num2str(t)],'Location','EastOutside')
function dzdt = odefun(t,z,y,sigma,beta,gamma)
xmy = abs(y-y');
dzdt = beta*(1-z).*trapz(y,exp(-sigma*xmy).*z',2)-gamma*z;
end
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!