Help with coding a model

1 visualizzazione (ultimi 30 giorni)
iheb kharab
iheb kharab il 9 Ott 2022
Modificato: iheb kharab il 25 Nov 2022
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

Risposta accettata

Davide Masiello
Davide Masiello il 9 Ott 2022
Modificato: Davide Masiello il 9 Ott 2022
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 Commento
iheb kharab
iheb kharab il 9 Ott 2022
Modificato: iheb kharab il 9 Ott 2022
Hi Davide,
Yes, this is exactly what I was looking for! Thank you so much.
Kind regards,
Iheb

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Mathematics in Help Center e File Exchange

Prodotti


Release

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by