Solving a system of differential equations with a variable stored in an array

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.

Answer 1

[Star Strider]

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)'])

.

[Mahima]

Whydo all 3 of the curves take a particular beta value for different time steps here? I was instead looking for them taking it at the same timestep. For e.g. we can clearly see that the color of beta hence the value is different for the last time step for all 3 of them.

[Sam Chak]

@Mahima, considering that there are three states, S, I, and R, are you suggesting to use only three distinct colors and assign 12 monochromatic color palettes to each color?

[Star Strider]

‘Whydo all 3 of the curves take a particular beta value for different time steps here?’

It is my understanding that is what you want them to do. I made you code a bit more efficient.

[Mahima]

@Star Strider what I mean is, let's take last 30 days for an example. I want the beta value to be same for all the 3 plots for those last 30 days. And I believe in your plot, same color means same value (sometimes it doesn't) but in your plot say in the last time step (last 30 days) the color is not same. Do those 3 plots have different values of beta for the same time step?

[Star Strider]

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)'])

.

Answer 2

[Torsten]

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)

[John D'Errico]

But this does not do as asked. (It may be more appropriate than what was asked, but that may not be pertinent.)

You use a linear interpolation on the vector b. However, it was requested that b change suddenly, every 30 days.

"I want the beta to change values every 30 days"

[Torsten]

Ok, result is the same:

b_array = [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 = 0.2;
x0 = [0.99999 0.00001 0];
hold on
for i = 1:12
  b = b_array(i);
  g = @(t,x)[-b*x(1)*x(2);b*x(1)*x(2)-k*x(2);k*x(2)];
  tspan = [(i-1)*30,i*30];
  [t,x] = ode45(g,tspan,x0);
  plot(t,x(:,1),'b',t,x(:,2),'r',t,x(:,3),'g')
  x0 = x(end,:);
end
hold off
grid on

[Mahima]

So, to get those 2 peaks what changes can help? I believe the paper changes k values also. Do you think iterating k values would help me get 2 clear peaks?

[Torsten]

I wonder what N is in your equations.

[Sam Chak]

Ha... I also overlooked this. 😅

Answer 3

[Sam Chak]

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)".

Edit: The code is updated to include N, as shown in the SIR model.

%% 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