This page contains a maple code for the course: MATH 111 Applied Mathematical Modelling I --- Maple Assignment V.
# endemic.maple Maple code to integrate a simple SIR model for an
# 14.10.03 endemic disease with vaccination.
# 10.10.07 Revised to run on maple 10.
# 24.12.15 Minor revisions to emphasis the final output point.
#
with(linalg):
with(plots):
R :=10.0; # the basic reproductive number.
vac:= 0.00; # the fraction of the population that is vaccinated.
tstart := 0; # the initial value for time.
tend := 10; # the final value for time.
# define the differential equations
de1 := diff(u(t),t) = b/(b+c)*(1-vac-u(t)) -R*u(t)*v(t);
de2 := diff(v(t),t) = (R*u(t)-1)*v(t);
de3 := diff(w(t),t) = b*vac/(b+c) +c/(c+b)*v(t) -b/(b+c)*w(t);
interval := 5.0: # number of points calculated per unit of time.
susceptibles := 0.98; # initial fraction of population that are susceptibles.
ic1 := u(0) = susceptibles; # initial pop'n fraction that are susceptibles.
ic2 := v(0) = 1-susceptibles; # initial pop'n fraction that are infected.
ic3 := w(0) = 0; # initial pop'n fraction that are immune.
deqs := {de1,de2,de3}: # define the system of DEs.
ics := {ic1,ic2,ic3}: # define the initial conditions.
vars := {u(t),v(t),w(t)}: # define the dependent variables.
b := 10.0; # per-capita birth rate.
c := 10*b; # the rate at which invected people recover
# from the disease.
# define the values of time at which the output is calculated.
ans := array([seq(i/interval,i=tstart..(interval*tend))]):
# integrate the system of differential equations
#
solnum := dsolve(deqs union ics,vars, \
type=numeric, output=ans, maxfun=300000):
# Print the values of the independent and dependent variables at the final
# integration point (tend). NOTE that the output is NOT necessarily in the
# format u(t), v(t), w(t).
#
# The next command pulls out an array that contains the time and the
# solution points.
solmatrix := eval(solnum[2,1]):
rowsize := rowdim(solmatrix):
print("Print the output at the final integration point - these numbers may be important!");
eval(solnum[1,1]); # this gives us the output order
# Print the output at the final integration point: important!
row(solmatrix,rowsize);
# Plot some figures. You WILL need to play around with these
# commands to get the best possible figures for your report.
#
# Figure One. Plot the solution in the susceptible-invected plane.
odeplot(solnum,[u(t),v(t)],labels=["u(t)","v(t)"]);
# Figure Two. Plot the susceptible-time and invected-time data on the
# same figure. Make sure you know which one is which.
p1 := odeplot(solnum,[t,u(t)],colour=BLUE):
p2 := odeplot(solnum,[t,v(t)],colour=RED):
display({p1,p2},title="Susceptibles and Infectives",\
labels=["time","fractional population"]);
# Figure Three. Plot the fraction of infected people in the population.
# You may find it useful to change the `y'-view to
# clarrify the long-term behaviour: is it heading towards zero?
odeplot(solnum,[t,v(t)],colour=RED,labels=["time","v(t)"],\
title="Fractional infected population",\
view=[tstart..tend,0..0.1]);
# Figure Four. Plot the immune-time data.
odeplot(solnum,[t,w(t)],labels=["t","w(t)"]);
R := 'R':
b := 'b':
c := 'c':
deqs := 'deqs':
de1 := 'de1':
de2 := 'de2':
de3 := 'de3':
ics := 'ics':
ic1 := 'ic1':
ic2 := 'ic2':
ic3 := 'ic3':
interval := 'interval':
output := 'output':
p1 := 'p1':
p2 := 'p2':
rowsize := 'rowsize':
solmatrix := 'solmatrix':
solnum := 'solnum':
susceptibles := 'susceptibles':
tend := 'tend':
tstart := 'tstart':
vac := 'vac':
vars := 'vars':