Description

Contents
1 Introduction 1
2 The problems 1
1 Introduction
This note contains the list of problems for our lab session
2 The problems
Problem 1
1. Copy script rode mwe1 in work/l7p1.m and adapt it to solving the initial value problem
y0(t) = 1 − y(t)2, y(0) = 0
for t ∈ [0,2].
2. For each value of t examine Richardson’s fraction. Are they behaving in a manner which is consistint with the existence of an asymptotic error expansion?
3. Show that y(t) = tanh(t) is the solution of this initial value problem and include this information in your script.
4. Is there good agreement between the error estimates and the true error?
Remark 1 The differential equation
y0(t) = 1 − y(t)2 (1)
(2)
g √
provided a = pk and b = gk. This is the differential equation which corresponds to a rock falling straight down through a homogeneous atmosphere.
Problem 2 Implementing and testing new method
1. Study the implementation of rk.m and at least one of the dependencies, say, phi4.m to the point were you can implement the method
) (3)
(4)
, (5)
(6)
as work/psi.m.
2. Develop a script work/l7p2.m which demonstrates that the global error for this method is O(h3).
3. Compare the error estimates to the exact error for a differential equation where you know the solution. Can you trust the error estimates?
Problem 3 The connection between regular integration and solving ordinary differential equations Consider the problem of computing the standard normal distribution function, i.e
(7)
1. Show that t → F(t) is the unique solution of the initial value problem
, (8)
(9)
where f : R×R→R given by
, (10)
is independent of y.
2. Develop a script work/l7p3.m which computes a table of t → F(t) for 21 equidistant points in the interval [0,2]. Your relative error must be less that τ = 10−7.
3. This goal can be accomplished using a smallest stepsize of h = 1/640. Which stepsize did you use?
Remark 2 At this point you can compute reliable error estimates. Regardless, the true solution is available in MATLAB as sol=@(t)0.5*(1+erf(t/sqrt(2))), where erf is MATLAB’s implementation of the standard error function.
Problem 4 The SIR model of infectious diseases. Infectious diseases such Ebola in can modelled using differential equations. In this problem we consider the standard SIR model for a epidemic from which everybody eventually recovers. It has three groups of people, susceptible (S), infected (I), and recovered (R). Susceptible people become infected at a rate which is proportional to the product of their number and the number of infected. Infected people recover at a constant rate and develop natural immunity. The corresponding system of ordinary differential equations is
S0 = −αIS, (11)
I0 = αIS − βI, (12)
R0 = βI, (13)
where α,β > 0 are constants which determine the rate of infection/rate of recovery.
1. Implement a function work/viral.m similar to shell4.m which models the above system of ordinary differential equations. The call sequence must be y=viral(t,x) where
(a) t is a dummy variable which must be present in order for viral.m to be compatible with the main program rk.m.
(b) x is a column vector, such that x(1) is the fraction of the population which is susceptible, x(2) is the fraction of the population which is infected, x(3) is the fraction of the population which has recovered and is immune.
(c) y is a column vector such that y(i) is the time derivative of x(i).
2. Develop a script work/l7p4.m which simulates an epidemic for which α = 0.5 days−1 and β = 0.04 days−1. Initially, 99% of the population is susceptible and 1% is infected. Track the disease using rk.m for 60 days. Verify, that more than 85% of the population has recovered by the end of the period.
3. An outbreak is said to be contained if the number of sick people is not increasing. Show that any outbreak will be contained if and only if .
4. Return to example and verify that the number of infected peaked at the time T where . Verify that more than 70% of the population was ill at the time!
5. The purpose of vaccination is ensure that outbreaks are always contained. If a fraction z of the population is immune to the disease, then the model changes to
S0 = −αI(1 − z)S (14)
I0 = αI(1 − z)S − βI (15)
R0 = βI (16)
Change your implementation of viral to contain a variable z which specifies the fraction of the population which is immune to the disease.
6. (Herd immunity) Show that outbreaks are always contained if and only if S(0) <
β
.
7. Return to the previous example, but assume that 95% of the population has be vaccinated at birth and track the outbreak for 120 days. Verify, that the disease can never threatend the fabric of society, because the sick are always so few in number that adequate care can be provide while normal functions continue.