Euler's Method

11 次查看（过去 30 天）
Edinson Manga
Edinson Manga 2020-6-7
评论： darova 2020-6-12
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).
  1 个评论
Edinson Manga
Edinson Manga 2020-6-7
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Alan Stevens
Alan Stevens 2020-6-7
编辑：Alan Stevens 2020-6-7
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)])
  4 个评论
显示 2更早的评论
Alan Stevens
Alan Stevens 2020-6-8
编辑：Alan Stevens 2020-6-8
What error? It works perfectly for me! Describe the error message that you get.
darova
darova 2020-6-12
i changed this line
Average = round(mean(y));
and it works ok for me too

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 App Building 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by