Euler's Method

Hi!

Model of an epidemic: We will now describe a mathematical model for the spread of an epidemic.

Suppose we have a community of L people that initially contains P infected people and Q uninfected people. Let y (t) be the number of people infected at one time t. If the disease is not very serious, like the common cold, everyone continues to be active and the epidemic spreads.

Since there are P Q possible contacts between people from one group to another, the rate of change of y (t) is proportional to

P Q, so the problem can be modeled using the initial value problem

(a) Taking L = 25000, k = 0.00003 and h = 0.2, with the initial condition y0 = 250, use

Euler's method to approximate the solution of the P.V.I. in the interval [0, 60].

(b) Draw the graph of the approximate solution in part (a).

(c) Estimate the average number of people infected by calculating the arithmetic mean

of section (a).

[Edinson Manga]

H = [0.2]; 
Y_euler = cell(1,numel(H));
for k=1:numel(H)
    h = H(k);
    x = 0:h:60;  
    y = zeros(1, numel(x)); 
    y(1) = 250;  
    n = numel(y); 
    for i=1:n-1
        f = ((0.00003*y(i))*(25000-y(i)));
        y(i+1) = y(i) + h * f;
    end
    Y_euler{k} = [x; y].';
end
syms Y(X)
eq = diff(Y) == (0.00003*Y)*(25000-Y);
ic = Y(0) == 250;
sol = dsolve(eq, ic);
y_sol = double(subs(sol, X, x));
hold on
plot(x, y_sol, 'r', 'DisplayName', 'true solution');
for k=1:numel(H)
    plot(Y_euler{k}(:,1), Y_euler{k}(:,2), 'b--', 'DisplayName', ['Euler method h=' num2str(H(k), '%.3f')]);
end
legend

This is my attempt, but I still can't find a way to find the arithmetic mean.

Answer 1

[Alan Stevens]

Here's a rather simpler way to do what you want:

% EulerEpidemic.m

% Data
L = 25000; k = 0.00003;  h = 0.2;
y0 = 250;  tend = 60;
n = tend/h;

% Initial conditions 
y(1) = y0;
t(1) = 0;

% Euler
for i = 2:n+1
    f = k*y(i-1)*(L - y(i-1));
    y(i) = y(i-1) + h*f;
    t(i) = (i-1)*h;
end

Average = round(mean(y),0);

% True
i = 1:n+1;
y_true(i) = L./( 1 - (1-L/y0)*exp(-k*L.*t(i)) );

% Graph
plot(t,y,'--',t,y_true,'r'),grid
xlabel('time'),ylabel('infected')
legend('Euler','True')
text(20,L/2,['Average = ', num2str(Average)])

[Alan Stevens]

Edited code to correct small error.

[Edinson Manga]

When I execute the code, I get an error

[Alan Stevens]

What error? It works perfectly for me! Describe the error message that you get.

[darova]

i changed this line

Average = round(mean(y));

and it works ok for me too