Help with coding a model

1 回表示 (過去 30 日間)
iheb kharab
iheb kharab 2022 年 10 月 9 日
編集済み: iheb kharab 2022 年 11 月 25 日
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

採用された回答

Davide Masiello
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
  1 件のコメント
iheb kharab
iheb kharab 2022 年 10 月 9 日
編集済み: iheb kharab 2022 年 10 月 9 日
Hi Davide,
Yes, this is exactly what I was looking for! Thank you so much.
Kind regards,
Iheb

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeMathematics についてさらに検索

製品


リリース

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by