Clear Filters
Clear Filters

Help with coding a model

1 view (last 30 days)
iheb kharab
iheb kharab on 9 Oct 2022
Edited: iheb kharab on 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

Accepted Answer

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

Sign in to comment.

More Answers (0)

Categories

Find more on Mathematics in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!