Solving a system of differential equations with a variable stored in an array
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
I have these 3 equations:
dS/dt = -β * S * I/N,
dI/dt = β * S * I/N - γ * I,
dR/dt = γ * I,
In my case, the beta is varying and I have an array of the values that beta is supposed to take
beta_values = [0.50, 0.49, 0.45, 0.33, 0.27, 0.32, 0.37, 0.47, 0.46, 0.50, 0.50, 0.50]
Now, I want to run the code for 360 days and I want the beta to change values every 30 days, hence 30*12 = 360.
So far, I've written this :
syms t x
b = [0.49*ones(1,30), 0.39*ones(1,30), 0.48*ones(1,30), 0.48*ones(1,30), 0.43*ones(1,30), 0.28*ones(1,30), 0.37*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30)];
k = 1/5
for j=1:360
g = @(t,x)[-b(j)*x(1)*x(2);b(j)*x(1)*x(2)-k*x(2);k*x(2)]
[t,xa] = ode45(@(t,x) g(t,x),[0 360],[0.99999 0.00001 0]);
end
figure
for i=1:3
plot(t,xa(:,i))
hold on;
title(['y(t) for a=',num2str(i)'])
end
This is giving me a plot, but that has just 1 peak and is not quite what I was looking for.
And this is what I want:
I'm new to MATLAB therefore not even sure where I went wrong.
Looking for help. Thanks.
0 comentarios
Respuestas (3)
Star Strider
el 19 de Mzo. de 2024
The β value does not cause the curves to vary much. There may be other problems with your initial conditions for the outbreaks or the way you set this up.
This varies β a bit more efficiently, and shows the results —
% syms t x
% b = [0.49*ones(1,30), 0.39*ones(1,30), 0.48*ones(1,30), 0.48*ones(1,30), 0.43*ones(1,30), 0.28*ones(1,30), 0.37*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30)];
beta = [0.50, 0.49, 0.45, 0.33, 0.27, 0.32, 0.37, 0.47, 0.46, 0.50, 0.50, 0.50];
k = 1/5;
g = @(t,x,b)[-b*x(1)*x(2);b*x(1)*x(2)-k*x(2);k*x(2)];
tspan = linspace(1, 30, 30);
ic = [0.99999 0.00001 0];
for j=1:12
[t{j},xa{j}] = ode45(@(t,x) g(t,x,beta(j)),[tspan+30*(j-1)], ic);
ic = xa{j}(end,:);
end
figure
for i=1:numel(t)
hp{i} = plot(t{i},xa{i}, 'DisplayName',sprintf('\\beta = %.2f',beta(i)));
hp1{i} = hp{i}(1);
hold on
end
hold off
grid
axis('padded')
xlabel('Time')
ylabel('Number of Cases')
legend([hp1{:}], 'Location','best')
% title(['y(t) for a=',num2str(i)'])
.
5 comentarios
Star Strider
el 20 de Mzo. de 2024
The colours in the curves have nothing to do with the β value, they relate to the order in which the curves are drawn. The β in a specific time (month) is the same for all the curves in that month.
% syms t x
% b = [0.49*ones(1,30), 0.39*ones(1,30), 0.48*ones(1,30), 0.48*ones(1,30), 0.43*ones(1,30), 0.28*ones(1,30), 0.37*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30), 0.50*ones(1,30)];
cm = colormap(turbo(12));
beta = [0.50, 0.49, 0.45, 0.33, 0.27, 0.32, 0.37, 0.47, 0.46, 0.50, 0.50, 0.50];
k = 1/5;
g = @(t,x,b)[-b*x(1)*x(2);b*x(1)*x(2)-k*x(2);k*x(2)];
tspan = linspace(1, 30, 30);
ic = [0.99999 0.00001 0];
for j=1:12
[t{j},xa{j}] = ode45(@(t,x) g(t,x,beta(j)),[tspan+30*(j-1)], ic);
ic = xa{j}(end,:);
end
figure
for i=1:numel(t)
hp{i} = plot(t{i},xa{i}, 'DisplayName',sprintf('\\beta = %.2f',beta(i)));
hp1{i} = hp{i}(1);
hold on
end
hold off
% grid
xline(0:30:360, ':k', 'LineWidth',1.3)
axis('padded')
xlabel('Time')
ylabel('Number of Cases')
legend([hp1{:}], 'Location','best')
% title(['y(t) for a=',num2str(i)'])
.
Torsten
el 19 de Mzo. de 2024
This is what you get with your equations and your data for beta:
T = 0:30:360;
b = [0.50, 0.49, 0.45, 0.33, 0.27, 0.32, 0.37, 0.47, 0.46, 0.50, 0.50, 0.50];
b = [b(1) (b(1:end-1)+b(2:end))/2 b(end)];
k = 1/5;
B = @(t)interp1(T,b,t);
g = @(t,x)[-B(t)*x(1)*x(2);B(t)*x(1)*x(2)-k*x(2);k*x(2)];
[t,xa] = ode45(@(t,x) g(t,x),[0 360],[0.99999 0.00001 0]);
plot(t,xa)
5 comentarios
Sam Chak
el 19 de Mzo. de 2024
Editada: Sam Chak
el 20 de Mzo. de 2024
Hi @Mahima
In addition to the standard integration methods suggested by @Star Strider and @Torsten, you can solve the SIR model by mathematically describing the piecewise beta function. It's akin to stepping stones reminiscent of the Giant's Causeway in Ireland. You can use a formula that I fondly refer to as the "Piecewise Function Put Together (PFPT)".
%% Call ode45
T = 0:30:360;
t = linspace(0, 360, 36001);
[t, xa] = ode45(@sirODE, t, [0.99999 0.00001 0]);
B = zeros(1, numel(t));
for j = 1:numel(t)
[~, B(j)] = sirODE(t(j), xa(j,:).');
end
%% Plot results
subplot(211)
plot(t, B, 'linewidth', 1.5, 'Color', '#265EF5'), grid on, axis([0, 360, 0, 0.6])
xticks(T), ylabel \beta, title('\beta')
subplot(212)
plot(t, xa), grid on, xlim([0, 360])
xticks(T), xlabel Days, title('SIR'), legend('S', 'I', 'R', 'location', 'east')
%% SIR Model with the piecewise beta function
function [dx, B] = sirODE(t, x)
S = x(1);
I = x(2);
R = x(3);
N = x(1) + x(2) + x(3);
b = [0.50, 0.49, 0.45, 0.33, 0.27, 0.32, 0.37, 0.47, 0.46, 0.50, 0.50, 0.50];
k = 1/5;
T = 0:30:360;
M = zeros(numel(b), numel(t));
for i = 1:numel(b)
M(i,:) = heaviside(t - T(i));
end
% Construction of the beta function using PFPT formula
B = b(1)*M(1,:) - (b(1) - b(2))*M(2,:) - (b(2) - b(3))*M(3,:) - (b(3) - b(4))*M(4,:) - (b(4) - b(5))*M(5,:) - (b(5) - b(6))*M(6,:) - (b(6) - b(7))*M(7,:) - (b(7) - b(8))*M(8,:) - (b(8) - b(9))*M(9,:) - (b(9) - b(10))*M(10,:) - (b(10) - b(11))*M(11,:) - (b(11) - b(12))*M(12,:);
% ODEs of SIR model
dx = [-B*S*I/N;
B*S*I/N - k*I;
k*I];
end
0 comentarios
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!