I need to make a code to get phase portraits for a mathematical model with "history". I will explain after the code.
close all;
clear all;
times = 1990:1:2015;
hold on
b=zeros(1,26); %75-2000 per 5 years
b(1:5)=0.0358;
b(6:10)=0.0339;
b(11:15)=0.0311;
b(16:20)=0.0275;
b(21:26)=0.0249;
m=zeros(1,26); %90-2015 per 5 years
m(1:5)=0.008;
m(6:10)=0.0031;
m(11:15)=0.0137;
m(16:20)=0.0147;
m(21:26)=0.0125;
l=zeros(1,26); %90-2015 per 5 years
l(1:5)=0.015;
l(6:10)=0.031;
l(11:15)=0.026;
l(16:20)=0.015;
l(21:26)=0.014;
u=zeros(1,26); %90-2015 per 5 years
u(1:5)=0.04;
u(6:10)=0.02;
u(11:15)=0.038;
u(16:20)=0.05;
u(21:26)=0.035;
S=zeros(1,26);
I=zeros(1,26);
N=zeros(1,26);
S(1)=18442000;
I(1)=186000; %1990
N(1)=18628000;
P=zeros(1,26); %15 years before S
P(1:5)=12788000;
P(6:10)=14731000;
P(11:15)=16968000;
P(16:20)=19696000;
P(21:26)=22893000;
for i=1:26
[time, xy] = ode45('test_func',times,[S(i) I(i) N(i) P(i) b(i) m(i) l(i) u(i)]);
plot(time,xy(:,1),'-g',time,xy(:,2),'-r',time,xy(:,3),'-b');
end
function rhs = test_func(t,xx)
S = xx(1);
I = xx(2);
N = xx(3);
P = xx(4);
b = xx(5);
m = xx(6);
l = xx(7);
u = xx(8);
Sdot=b*P-m*S-l*S*I;
Idot=l*S*I-(m+u)*I;
Ndot=Sdot+Idot;
rhs = [Sdot; Idot; Ndot; P; b; m; l; u;];
end
List of details:
S=healthy populationI=infected populationN=total populationb=birth rate (15 years prior to S)m=mortality ratel=infection chance on contactu=death rate due to diseaseP and S represent the same thing just in different time periods (P=15 years prior to S), also all P values are given.
The code needs to return the phase portrait of S,I and N. I am definitely not 100% sure my code is right for what I aim to do but this is what I came up with. Currently the code runs but never ends. Any suggestions on my code or help to fix the error are welcome.
I was also thinking of adding the following inside the for loop right between the ode45 and plot, if necessary:
if i<26
xy(i+1)=S(i+1);
xy(i+27)=I(i+1);
xy(i+53)=N(i+1);
end
ode45 is intended to solve ordinary differential equations. An ordinary differential equation problem setup will have some essential components.
xdot = f(t, x)
is the differential equation.
x(t=t0) = x0
is the initial condition.
t is the independent variable and corresponds to time in your implementation.
x is the dependent variable and corresponds to S, I, and N in your implentation.
The t0 and x0 in the initial condition correspond with times(1)=1990 and S(1), I(1), and N(1).
The remaining task is to define f in a way that MATLAB understands. Once you have all these components, ode45 is ready to use.
Defining f is probably the most difficult part. In your implementation, it corresponds with test_func and the additional parameters required for test_func (P, b, m, l, u).
It's important to note that in your situation, these parameters also depend on time. Perhaps it's clearer to write them as P(t), b(t), m(t), l(t), and u(t).
In this case, it's probably helpful to take a look at the interp1 function, which is a built-in linear interpolation function. Given your data points, MATLAB can estimate the value of P(t), b(t), m(t), l(t), and u(t), when t isn't every 5th year.
function q41895153_ode45()
time_range = 1990:0.1:2015;
time_hist = 1990:5:2015;
b=zeros(1,6); %75-2000 per 5 years
b(1)=0.0358;
b(2)=0.0339;
b(3)=0.0311;
b(4)=0.0275;
b(5)=0.0249;
b(6)=0.0249;
m=zeros(1,6); %90-2015 per 5 years
m(1)=0.008;
m(2)=0.0031;
m(3)=0.0137;
m(4)=0.0147;
m(5)=0.0125;
m(6)=0.0125;
l=zeros(1,6); %90-2015 per 5 years
%{
l(1)=0.015;
l(2)=0.031;
l(3)=0.026;
l(4)=0.015;
l(5)=0.014;
l(6)=0.014;
%}
l(1)=0.001;
l(2)=0.001;
l(3)=0.001;
l(4)=0.001;
l(5)=0.001;
l(6)=0.001;
u=zeros(1,6); %90-2015 per 5 years
u(1)=0.04;
u(2)=0.02;
u(3)=0.038;
u(4)=0.05;
u(5)=0.035;
u(6)=0.035;
S0=18442;
I0=186; %1990
P=zeros(1,6); %15 years before S
P(1)=12788;
P(2)=14731;
P(3)=16968;
P(4)=19696;
P(5)=22893;
P(6)=22893;
[time, xy] = ode45(@test_func,time_range,[S0 I0],odeset(),time_hist,P,b,m,l,u);
S = xy(:,1)
I = xy(:,2)
N = S + I
plot(time,xy);
end
function rhs = test_func(t,xx,time_hist,P,b,m,l,u)
S = xx(1);
I = xx(2);
% Interpolate to find b(t), m(t), l(t), u(t), P(t)
bt = interp1(time_hist,b,t);
mt = interp1(time_hist,m,t);
lt = interp1(time_hist,l,t);
ut = interp1(time_hist,u,t);
Pt = interp1(time_hist,P,t);
Sdot=bt*Pt-mt*S-lt*S*I;
Idot=lt*S*I-(mt+ut)*I;
rhs = [Sdot; Idot];
end