MA50174 ADVANCED NUMERICAL

METHODS
[Link]
Contents
52
Chapter 6
Initial Value Problems (IVPs)
6.1 Introduction
An initial value problem is an ordinary dierential equation of the form
du
dt
= f(t, u) where u(t
o
) = u
o
is given and u IR
n
. (6.1)
In the special case when f (t, u) = f(t) we have
u(t) =
_
t
to
f dt +u
o
,
and the task of evaluating this integral accurately is called quadrature To solve any dierential equation
we need to put it into the standard form given by (6.1). Any equation involving higher derivatives
can be reformulated as such a vector equation. For example, if w satises the second order dierential
equation,
d
2
w
dt
2
= w
Let
u
1
= w, u
2
= dw/dt.
Then
du
1
/dt = u
2
,
du
2
/dt = u
1
.
Thus if
u = (u
1
, u
2
)
T
we have
du
dt
=
_
0 1
1 0
_
u Au
Initial value problems (IVPs) come in various forms and there is no such thing as a perfect all purpose
IVP solver. MATLAB oers you quite a choice. Will try to show you how to choose which one to use
for a given problem.
We say that an ODE problem is
1. Linear if f (t, u) is linear in u
2. Autonomous if f (t, u) f (u)
3. Non-sti if all components of the equation evolve on the same timescale. This occurs (roughly)
if the Jacobian matrix f /u has all its eigenvalues of similar size.
53
4. Sti if dierent components of the system evolve on dierent time scales. These are very common
in chemical reactions with reactions going on at dierent rates and in ODEs resulting from spatial
discretisations of PDEs. They also occur in PDEs where dierent modes (Fourier modes) evolve
at very dierent rates. Sti problems are much harder to solve numerically than non-sti ones.
5. Hamiltonian if f takes the form
f = J
1
H
where H(y) is the Hamiltonian of the system and
J =
_
0 I
I 0
_
where I is the
n
2

n
2
identity matrix.
Hamiltonian equations arise very commonly in rigid body mechanics, celestial mechanics (astron-
omy) and molecular dynamics. To solve a Hamiltonian equation accurately over long time periods
we must use special numerical methods, such as symplectic or reversible methods.
All modern software for IVPs is a combination of three components
The actual solver.
A way of estimating the error of the solution.
A step-size control mechanism.
Thus, the IVP solver attempts to use the best method to keep the (estimated) error within a prescribed
tolerance. Whilst it is essential to have some form of error control, no such method is infallible. The
numerical solution of IVPs is well covered in many texts, for example [Link] A rst course in the
numerical solution of dierential equations.
6.2 Quadrature
Suppose that f(t) is an arbitrary function, how accurately can we nd
u =
_
b
a
f(t)dt ? (6.2)
MATLAB determines this integral approximately by using a composite Simpsons [Link] idea behind
this is as follows:
Suppose we take an interval of length 2h, without loss of generality this is the interval [0, 2h].
Evaluate f at the points 0, h, 2h
Then approximate f by a parabola through these points, by using a quadratic interpolant as in
Chapter 3.
Integrate this approximation to get an estimate for the integral of f.
This gives
_
2h
0
fdx
h
6
[f(0) + 4f(h) +f(2h)] S
h
This approximation is unreasonably accurate. It can be shown (see Froberg, Numerical Analysis)
that
|S
h

_
2h
0
f(t)dt| =
h
5
90
|f
(iv)
()| where 0 < < 2h (6.3)
54
The local error is proportional to h
5
and to f
(iv)
. The traditional use of Simpsons rule to evaluate (6.2)
over [a, b] breaks this interval into sub-intervals of length 2h takes h constant between a and b and adds
up the results to get a total error estimated by
h
4
90
|b a| max(|f
(iv)
|).
MATLAB is more intelligent than this and it uses an adaptive version of Simpsons rule. In this
procedure h is chosen carefully over each interval to keep the error estimated by (6.2) less than a user
specied tolerance. In particular h is small when f is varying more rapidly and f
(iv)
is large. The
integrals over each sub-interval are then combined to give the total. This method is especially eective
if f has a singularity.
The procedure uses the instruction
> i = quad(@fun, a, b, tol)
where fun is the function to be integrated and tol is the tolerance. There is another MATLAB code
> quadl. This approximates f by a higher order polynomial. This is more accurate if f is smooth but
unreliable if f has singularities.
6.3 Non-sti ordinary dierential equations (ODEs)
A non-sti ordinary dierential equation has components which all evolve on similar time-scales. They
are the easiest dierential equations to solve by using a numerical method. In particular they can often
be solved by using explicit methods that do not require the solution of nonlinear equations.
6.3.1 The Forward Euler Method
The oldest, easiest to apply and analyse, method for such problems is the explicit forward Euler method.
Suppose that u satises the ODE
du
dt
= f(t, u), u(0) = u
0
We take a small step h and approximate u((n 1)h) by U
n
We set U
1
= u
0
Now for each successive n we update U
n
through
U
n+1
= U
n
+hf (t, U
n
) = U
n
+hf ((n 1)h, U
n
) (6.4)
This method is easy to use and each step is fast as no equations need to be evaluated and there is only
one function evaluation per step. The Forward Euler method is still used when f is hard to evaluate and
there are a large number of simultaneous equations. Problems of this kind arise in weather forecasting.
The main problem with this method is that there are often severe restrictions on the size of h.
Local Error.
At each stage of the method a small error is made. These errors accumulate over successive intervals
to give an overall error. If the exact solution u(t) is substituted into the equation (6.4), there is a
mismatch E between the two sides, called the local truncation error or LTE. This is a good estimate
for the error made by the method at each step. It is estimated by
|E| = |u(nh)u((n1)h)hf((n1)h, u((n1)h))| =
h
2
2
|u

()| where (n1)h < < nh (6.5)

As in quadrature the local truncation error is proportional to a power of h, in this case h
2
and a higher
derivative of u, in this case u

. So, if u

is large (rapid change), we must take h small to give a small

55
error. The smaller the value of h is, the more accurate the answer will be. (We will see later that h is
also restricted by stability considerations).
Global Error
Each time the method is applied an error is made and the errors accumulate over all the calculations.
If we want to approximate u(T) then need to make T/h + 1 N calculations so that U
N
u(T). We
will call the global error if
= |U
N
u(T)|, N =
T
h
+ 1.
This can be crudely estimated by

Nh
2
2
max|u

| Th max|u

|.
We note that the overall error is proportional to h and to max|u

|. We say that this is an O(h) or a

rst order method.
The estimate we have obtained implies that the global error grows in proportion to both h and T. In
fact, this is rarely observed. In the worst case for some problems the error grows in proportion to e
T
,
when looking at periodic problems it can grow like T
2
and for problems with a lot of symmetry, the
errors can accumulate much more slowly so that the overall error may not grow at all with T.
We start by looking at the worst case analysis.
We say that f is globally Lipshitz if there is a constant L such that
|f(t, y) f(t, z)| < L|y z| (6.6)
It can then be shown (see A. Iserles A rst course in the numerical solution of dierential equations)
that if |u

| M then the global error has an upper bound given by

hM
2L
_
e
LT
1
_
(6.7)
If T is xed and h 0 then the upper bound given in (6.7) shows that is proportional to h in this
limit. In practice the Forward Euler method is rarely used in this way. We often take h xed and let
T . This is very dangerous as errors can accumulate rapidly making the method unusable. In
this case we need a more careful control on the global growth of errors. This is the philosophy behind
geometric integration based methods.
Variable time-stepping
As in the quadrature process, Eulers method can be used with a variable step size. This allows some
control over the local error. To do this we take a sequence of (small) time steps, so that the solution is
approximated at times t
k
with so that U
k
u(t
k
)
h
k
t
k+1
t
k
.
If we have an estimate for u

(), [t
k
, t
k+1
] then h
k
is chosen so that at each step
E =
h
2
k
|u

|
2
< Tol
where TOL is specied by the user. Of course, without knowing u(t) we do not know u

in general.
We can either estimate an upper bound a -priori to the calculation or we can try to deduce its value
(a-posteriori) from the numerical calculation. This latter procedure often uses an error estimate call
the Milne device. (See Iserles.)
56
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Example 1: Use of the Forward Euler Method %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Consider the ODE du/dt = -2t u^2, u(0) = 1 %
% ========================== %
% %
% This has the solution u(t) = 1/(1 + t^2) %
% ================== %
% %
% So that u(1) = 1/2 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% We can solve this using the Forward Euler method with %
% step size h, to find a numerical approximation U for u(1) %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Do a series of runs with h reducing in size
%
h = 1;
j = 1;
while j < 8
h = h/2;
%
% Number of time steps
%
N = 1/h + 1;
U = zeros(1,N);
t = zeros(1,N);
U(1) = 1;
t(1) = 0;
for i=2:N
%
% Euler step
%
57
t(i) = h + t(i-1);
U(i) = U(i-1) - 2*h*t(i-1)*U(i-1)^2;
end
hh(j) = h;
%
% Error
%
er(j) = U(N)-1/2;
j = j+1;
end
A = [hh er]
plot(hh,er)
-----------------------------------------------------------------------------------
%
% The resulting solution has the following form
%
>> eul
A =
0.5000 -0.2500
0.2500 -0.1209
0.1250 -0.0591
0.0625 -0.0293
0.0312 -0.0146
0.0156 -0.0073
0.0078 -0.0036
>> h Error
We can see from this that when h is halved, so is the error. The resulting solution for two values of h
is given below.
6.3.2 The Runge-Kutta method.
A much more locally accurate method is the Runge-Kutta method. Its most famous form is called the
explicit fourth order Runge-Kutta or the RK4 method. Suppose that the ODE is
du
dt
= f (t, u). Then if
58
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
h=0.25
h = 0.0078
Exact
t
u
Figure 6.1: Solution of the Forward Euler Method
0 0.1 0.2 0.3 0.4 0.5
0.25
0.2
0.15
0.1
0.05
0
h
Error
Figure 6.2: Error of the Forward Euler method as a function of h
we know U
n
, and set t = (n 1)h the value of U
n+1
is given by the sequence of operations
k
1
= hf (t, U
n
)
k
2
= hf
_
t +
h
2
, U
n
+
k
1
2
_
k
3
= hf
_
t +
h
2
, U
n
+
k
2
2
_
k
4
= hf (t +h, U
n
+k
3
)
U
n+1
= U
n
+
1
6
(k
1
+ 2k
2
+ 2k
3
+k
4
)
It can be shown that there is a value C which depends on f in a complex way such that a local truncation
error E = |u(t +h) U
n+1
| is bounded by
E Ch
5
Over a large number of steps these errors accumulate as before to give a global error of the form:
C(e
LT
1)h
4
The error is proportional to h
4
. Hence the name an order 4 method. The error for a given h is much
smaller than for the Forward Euler method. This method is VERY widely favoured as
59
1. It is easy to use and no equations need to be solved at each stage.
2. It is highly accurate for moderate h values
3. It is a one step method i.e. U
n+1
only depends on U
n
4. It is easy to start and easy to code.
For many people it is the ONLY method they ever use!! However, it has certain disadvantages.
1. The function f must be evaluated four times at each iteration. This may be dicult if f is hard
or expensive to evaluate.
2. Errors accumulate rapidly as T increases.
3. The method cannot be used for sti problems unless h is very small.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Example 2: Use of the 4th Order Runge-Kutta method %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Consider the ODE du/dt = -2t u^2, u(0) = 1 %
% ========================== %
% %
% This has the solution u(t) = 1/(1 + t^2) %
% ================== %
% %
% So that u(1) = 1/2 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% We can solve this using the RK4 method with %
% step size h, to find a numerical approximation for u(1) %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% The function is calculated in frk.m %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------------------------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Do a series of runs with h reducing in size %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60
h = 1;
j = 1;
while j < 8
h = h/2;
%
% Number of time steps
%
N = 1/h + 1;
U = zeros(1,N);
t = zeros(1,N);
%
% Initial values
%
u(1) = 1;
t(1) = 0;
for i=2:N
%
% Runge-Kutta step
%
th = h/2 + t(i-1);
t(i) = h + t(i-1);
k1 = h*frk(t(i-1),U(i-1));
k2 = h*frk(th,U(i-1)+k1/2);
k3 = h*frk(th,U(i-1)+k2/2);
k4 = h*frk(t(i),U(i-1)+k3);
U(i) = U(i-1) + (k1+2*k2+2*k3+k4)/6;
end
hh(j) = h;
%
% Error
%
er(j) = U(N)-1/2;
j = j+1;
end
61
rat = er./hh.^4;
A = [hh er rat]
plot(log(hh),log(abs(er)))
>> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>> %
>> % Test of the Runge-Kutta code
>> %
>> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>>
>> rung
A =
0.50000000000000 -0.00029847713504 -0.00477563416071
0.25000000000000 0.00001355253692 0.00346944945065
0.12500000000000 0.00000139255165 0.00570389155200
0.06250000000000 0.00000009811779 0.00643024720193
0.03125000000000 0.00000000640084 0.00671176833566
0.01562500000000 0.00000000040734 0.00683396495879
0.00781250000000 0.00000000002567 0.00689065456390
>> % h error error/h^4
>> %
>> % Note that error/h^4 is almost constant
>> %
>> diary off
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.5
0.6
0.7
0.8
0.9
1
t
u
I like to think of the RK4 method as being like a Ford Fiesta. It is easy to use, everyone uses it and it is
good value for money. If you are going to the shops (i.e. solving a relatively straight forward problem)
it is the method to use. However, it wont handle rough country nor would you want it for a very long
drive.
To improve the accuracy the RK4 method is often used together with another higher order method to
estimate the local error and choose h accordingly. Suppose that the value of U
n+1
is given by an RK4
method. We could also use a dierent higher order method to calculate a separate value

U
n+1
. Both
62
U
n+1
and

U
n+1
approximate u(t +h) so that:
U
n+1
= u(t +h) +Ch
5

U
n+1
= u(t +h) +Dh
6
Subtracting these estimates we have
U
n+1

U
n+1
|C|h
5
+|D|h
6
If h is small then |C|h
5
U
n+1

U
n+1
so the dierence between the two calculations gives an
estimate for the local error of the RK4 method.
We can now use this estimate as follows. Given U
n
and h
1. Using U
n
and step size h, calculate U
n+1
,

U
n+1
and E
n+1
= U
n+1

U
n+1
.
2. If
TOL
32
< E
n+1
< TOL then accept the step
3. If E
n+1
<
TOL
32
then set h = 2h and repeat from 1.
4. If E
n+1
> TOL then set h =
h
2
and repeat from 1.
This method keeps the local errors below the specied tolerance and also (due to 3) makes an ecient
choice of step size. However it does not control the growth of errors. MATLAB uses such a Runge-
Kutta 4,5 pair above the form developed by Dormand & Prince in the ode45 routine. There is a similar
(cheaper but less accurate) Runge-Kutta 2,3 pair implemented in the routine ode23. In this routine
you can set both the Absolute Tolerance (as above) or the Relative Tolerance TOL/ u . You can the
ode45 routine as follows:
>[t,U] = ode45 (@fun, trange, u_0, options)
Here t is the vector of times at which the approximate solution vector U is given. Because the step-size
h is chosen at each stage of the algorithm these times will not (necessarily) be equally spaced. The
function f(t, u) is specied by fun and u
0
is the vector of initial conditions. The time interval for the
calculation is given by trange. Setting
trange = [a b]
means that the ode45 routine will integrate the ODE between a and b, giving its output at the (variable)
time steps it computes. Alternatively
trange = [a: c: b]
leads to output at the points a, a+c, a+2c etc. The options command is optional, but it allows control
over the operation of the ode45 routine. In particular you can set the tolerances and/or request statistics
on the solution. The options are set by using the odeset routine. For example if you want an absolute
tolerance of 10
12
you type
>options = odeset (AbsTol, 1e^(-12}))
By default the absolute tolerance is 10
6
and the relative tolerance is 10
3
.
63
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Example 3: Use of ode45 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Consider the ODE du/dt = -2t u^2, u(0) = 1 %
% ========================== %
% %
% This has the solution u(t) = 1/(1 + t^2) %
% ================== %
% %
% So that u(1) = 1/2 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% We can solve this using the ode45 method with %
% tolerances 1e-12,to find a numerical approximation for u(1)%
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options=odeset(AbsTol,1e-12,RelTol,1e-12);
trange = [0:0.05:1];
unit = [1];
[t,U] = ode45(@frk,trange,unit,options);
s = size(U);
err = U(s(1)) - 0.5
6.3.3 The St ormer - Verlet method and geometric integration
Geometric integration is a branch of numerical analysis which aims in part to control global error growth
of a numerical approximation over long times. An important application of geometric methods is to
systems of the form
u = v
v = f(u)
_
_
u +f(u) = 0.
More generally, geometric methods can be applied to Hamiltonian systems for which
q = H/p, p = H/q
For the problem u +f(u) = 0 we set q = u, p = du/dt and
H = p
2
/2 +F(q) with F =
_
fdq
64
Multiplying the dierential equation by u and integrating with respect to time we nd that
H = u
2
/2 +F(u) = const.
More generally, in an autonomous Hamiltonian system, the Hamiltonian H is a constant for all times.
This is an example of a conservation law. Many physical systems conserve particular quantities over
all times. An excellent example of this is the solar system considered in isolation with the rest of the
universe. This obeys a complicated set of dierential equations with complex (and indeed chaotic)
solutions. However the total energy, angular momentum and linear momentum are conserved for all
time.
To retain the correct dynamics of such a system in a numerical approximation it is essential that the
numerical approximation either exactly conserves the same invariants or (more usually) the approximate
equivalent of any conserved quantity varies from a constant by a small but bounded amount for all times.
If this occurs then the approximate solution is likely to be much closer to the true solution for all times.
The Forward Euler method, RK4 and ode45 are not good at preserving such conserved quantities. For
example if we consider the system
u +f(u) = 0,
u
2
2
+F(u) = H
An application of the Forward Euler method with U
n
u((n 1)h), V
n
u((n 1)h) gives
U
n+1
= U
n
+hV
n
, V
n+1
= V
n
hf(U
n
)
so that
H
n+1
=
(V
n+1
)
2
2
+F(U
n+1
)
= [V
n
hf(U
n
)]
2
/2 +F(U
n
+hV
n
)
= H
n
+
h
2
2
_
f(U
n
)
2
+ (V
n
)
2
f/u
_
+O(h
3
)
Thus H changes by
H
n+1
H
n
=
h
2
2
_
f(U
n
)
2
+ (V
n
)
2
f/u
_
+O(h
3
)
at each iteration. If f/u > 0 then H is not conserved and increases at each iteration as the errors
accumulate so that the global error in H grows like nh
2
. In the RK4 method we have the better result
that
H
n+1
H
n
= h
5
C
n
.
Again, in general the values of C
n
are positive for many problems and the errors in H accumulate
over a large number of iterations, with global errors growing like nh
5
. A widely used method which
avoids many of these problems and has excellent conservation problems is the St ormer - Verlet method
(SV) . This is the method of choice for simulations of celestial and molecular dynamics. Not only is
it (much) better than RK4 for long-time integrations it is also much cheaper and easier to code up as
it only requires one function evaluation per time step. It is also explicit, has global errors O(h
2
) and
it is symmetric. The disadvantage of the SV method is that it can only be used for a certain class
of problems (those with a separable Hamiltonian) and it is not a black-box code i.e. it requires some
thought to use it. However this is precisely what mathematicians are paid to do!
For the problem u +f(u) = 0 the SV method takes the following form
U

= U
n
+
h
2
V
n
V
n+1
= V
n
hf(U

)
U
n+1
= U

+
h
2
V
n+1
In the SV method we have
H
n+1
H
n
= h
3
D
n
.
65
This appears to be worse than RK4. However, unlike RK4 the errors DO NOT accumulate and tend
to cancel out. It can be shown that if H is the exact Hamiltonian then
|H
n
H| < Dh
3
where D does not depend on n. So, although H is not exactly conserved, the method stays close to it
for all time.
We now apply this to an example with f(u) = u
3
.
66
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Example 4: Stormer-Verlet method for %
% %
% u + u^3 = 0, u(0) = 1, u(0) = 0 %
% %
% In the exact eqn H is constant where %
% %
% H = (u)^2/2 + u^4/4 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t(1) = 0;
U(1) = 1;
V(1) = 0;
h = 0.1;
H(1) = 1/4;
for i=2:100
Usta = U(i-1) + (h/2)*V(i-1);
V(i) = V(i-1) - h*Usta^3;
U(i) = Usta + (h/2)*V(i);
t(i) = t(i-1) + h;
H(i) = V(i)^2/2 + U(i)^4/4;
end
plot(t,H)
67
1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1
1
0.8
0.6
0.4
0.2
0
0.2
0.4
0.6
0.8
1
u
v
Figure 6.3: A plot of U, V for the exact and numerical solution
0 1 2 3 4 5 6 7 8 9 10
0.25
0.252
t
H
Figure 6.4: A plot of H showing the bounded variation
68
6.4 Sti Dierential Equations
6.4.1 Denition
In a sti dierential equation the solution components evolve on very dierent timescales. This causes
a problem as a numerical method, such as ode45, possibly chooses a step size for the most rapidly
evolving component, even if its contribution to the solution is negligable. This leads to very small step
sizes, highly inecient computations and long waits for the user! The reason for this is an instability
in the method , where a small error may grow rapidly with each step.
Suppose we want to solve the ODE
u = f (u)
and the numerical method makes a small error e. We ask the question, how does e grow during the
calculation? To answer this we start by looking at how a small disturbance to the solution of the ODE
changes. Suppose that e is such a disturbance so that
d
dt
(u +e) = f (u +e)
To leading order the perturbation e then satises the linear dierential equation
du
dt
+
de
dt
= f (u) +Ae A =
f
u
.
The perturbation growth is thus described by the dierential equation
de
dt
= Ae where A = f /u. (6.8)
Now look at the ODE (6.8). This has the solution
e = e
At
e
0
where e
0
is the initial perturbation. So that if A = UU
1
then
e = Ue
t
U
1
e
o
.
Alternatively, if A has eigenvalues
i
and eigenvectors
i
then
A
i
=
i

i

i
: eigenvector
Thus if e
0
is in the direction of
i
so that
e
0
= ai
it follows that
e(t) = ae

i
t

i
It follows further that
e(t) = |a|e

i
t

i

where
i
is the real part of
i
so that it is the real part of the eigenvalues which control the growth of
small perturbations.
If |
i
| is large and the real part of
i
is less than zero, then the contribution to e in the direction of
the eigenvector
i
rapidly decays to zero. After only a short time the perturbation is dominated by
the component in the direction of the eigenvector
j
for which
j
has the largest real part over all the
eigenvalues.
We are now able to dene what we mean by a sti system.
DEFINITION
The system is sti if the matrix A has eigenvalues
i
for which
max
j
|
j
| min
j
|
j
|.
69
Typically, in an application a ratio of over 10 is considered to lead to a sti system
In a physical system the components for which |
j
| is large, and the real part of
j
is negative, decay
rapidly and are not seen in the solution apart from some initial transients. It is therefore somewhat
paradoxical that it is precisely these components which lead to instabilities in the numerical [Link]
is another way of saying that the solution has components that evolve at very dierent rates.
6.4.2 Numerical Methods
Now we look at a variety of numerical methods for solving the linear equation (6.8) so that we may
compare the growth of perturbations to the numerical solution with those of the true solution. First
look at the performance of the Forward Euler method when applied to this problem. We have
U
n+1
= U
n
+hf(U
n
)
so that
U
n+1
= U
n
+hAU
n
= (I +hA) U
n
Therefore
U
n
= (I +hA)
n
U
1
But I +hA has the same eigenvectors
j
as A and has eigenvalues
1 +h
j
So, the contribution to U
n
in the direction of the eigenvector
j
grows as
(1 +h
j
)
n1
or more precisely as |1 + h
j
|
n1
. In particular the errors caused by truncation error or by rounding
error only decay if |1 +h
j
| < 1 for all
j
.
We now nd an extraordinary paradox. Suppose that
j
is real and that

j
= with 1
Then
e

j
t
= e
t
1, if t 1
However
|1 +h
j
|
n
= |1 h|
n
1 if h > 2 and n 1.
So the most rapidly decaying components of the perturbation to the continuous solution are the com-
ponents of U
n
which are growing most rapidly.
DEFINITION
We say the Forward Euler method with step size h is stable if given a matrix A with eigenvalues
j
with the real part less than zero then |1 + h
j
| < 1 for all
j
. In particular if
j
are all real then the
method is stable only if
h < 2/max|
j
|
This is the restriction on h which means that solution errors do not grow. If h is larger than this bound
then errors grow and the method is unstable.
Here we see the problem. A component in the direction of
j
with large |
j
| and with real part
j
< 0
dominates the choice of step size even though this component is very small.
70
The LOCAL (TRUNCATION) ERROR that is made at each stage of the calculation is given by
h
2
|u

|
2
.
So the error made at each stage depends on the SOLUTION. If all of the eigenvalues of A have
negative real part then this error is ultimately dominated by the component of the solution that
decays most slowly, which is in turn determinated by the eigenvalue of A with the largest real part.
Soi that it is the eigenvalue with the real part closest to zero, and typically this is the eigenvalue
with the smallest modulus.
In contrast the GROWTH of the error as the iteration proceeds depends on the eigenvalue of A
with the largest modulus, even though this eigenvalue may have a large negative real part and
thus not contribute to the solution in any way.
There are therefore two restrictions on h, it must be small both for accuracy at each stage and for
stability to stop the errors growing. Stiness arises when the restriction on h for stability is much more
severe than the restriction for accuracy.
This introduces us to the ideas of stability and instability. It is surprisingly hard to give a precise
denition of what we mean by instabiliy in a numerical method which accounts for all cases - later on we
will give a precise denition which covers certain cases, but for the present we will have the following.
INFORMAL DEFINITION OF INSTABILITY
A numerical method to solve a dierential equation is unstable if the errors it makes (or indeed its
solution) grow more rapidly than the underlying solution. If the errors decay then it is stable.
Exercise Check the denition of stability for the Forward Euler method is consistent with this informal
denition.
Returning to the Forward Euler example. If h > 0 and = p +iq then |1 +h| < 1 implies that
(1 +hp)
2
+ (hq)
2
< 1
2hp +h
2
p
2
+h
2
q
2
< 0
2p +hp
2
+hq
2
< 0.
So that the method is stable if (p, q) lies in the circle of radius
1
h
shown in Figure 6.5. In this gure
30 20 10 0 10 20 30
30
20
10
0
10
20
30
p
q
STABLE
UNSTABLE
Figure 6.5: Stability region for the Forward Euler method
the shaded region shows the values of p and q for which the numerical method is stable. Recall that
the original dierential equation is stable provided that p lies in the half-plane p < 0. The shaded
71
region only occupies a fraction of this half-plane, although the size of the shaded region increases as
h 0. Thus, for a xed value of h the numerical method will only have errors which do not grow if
the eigenvalues of A are severely constrained. Unfortunately this is often not the case, particularly in
discretisations of partial dierential equations.
6.4.3 The Backward-Euler Method... an A-stable sti solver
We now look at another method, The Backward Euler Method (also called the BDF1 method). This is
given by
U
n+1
= U
n
+hf (U
n+1
)
The backward Euler method is much harder to use than Forward Euler as we must solve an equation
(which is usually nonlinear) at each step of the calculation to nd U
n+1
. It has the same order of error.
i.e. the global error is proportional to h. If we now apply this to the equation u = Au we have
U
n+1
= U
n
+hAU
n+1
This is a linear system which we need to invert to give
U
n+1
= (I hA)
1
U
n
Now, if the eigenvalues of A are
j
, those of (I hA)
1
are (1 h
j
)
1
, with the same eigenvectors.
Exercise: Prove this.
Thus the contribution to U
n
in the direction of
j
evolves as
|1 h
j
|
1n
Now let
j
= p +iq as before. It follows that
|1 h
j
|
2
=
1
(1 hp)
2
+ (hq)
2
Therefore the errors decay and the method is stable if
1
(1 hp)
2
+ (hq)
2
< 1
if 1 < (1 hp)
2
+ (hq)
2
or 0 < h(p
2
+q
2
) 2p
The resulting stability region is illustrated (shaded) in Figure 6.6: This picture is in complete contrast
to the one that we obtained for the Forward Euler method. The stability region is now very large and
certainly includes the half-plane p < 0. Thus any errors in the numerical method will be rapidly damped
out. Unfortunately the numerical solution can decay even if p 0, so that neutral or growing terms
in the underlying solution can be damped out as well. This is a source of (potential) long term error,
especially in Hamiltonian problems.
The Backward Euler method is very reliable and has other nice properties (a maximum principle) which
make it especially suitable for solving PDEs. The main disadvantage to using it is that we must solve
a nonlinear equation to nd U
n+1
. This is an example of the basic principle that there is no such
thing as a free lunch! The penalty of a stable method is the need to do more work. Note, however,
that the St ormer - Verlet method is a good approximation to a free lunch. Finding U
n+1
when the
system has a high dimension (e.g.10
4
) is a considerable task, especially as this calculation must be done
quickly and often. Usually we use an iterative method such as the Newton-Raphson or Broyden method.
Fortunately a good initial guess is available namely a value

U
n+1
which is obtained by taking one-step
of an explicit method (e.g Forward Euler) applied to U
n
. This procedure is called a predictor-corrector
method in which the explicit method predicts the value

U
n+1
which is then corrected (usually by an
iterative method) to give U
n+1
.
72
30 20 10 0 10 20 30
30
20
10
0
10
20
30
p
q
STABLE
UNSTABLE
2/h
1/h
Figure 6.6: Stability region for the Backward Euler method.
6.4.4 The Trapezium Rule: a symmetric sti solver
The Trapezium Rule is given by
U
n+1
= U
n
+
h
2
[f (U
n
) +f (U
n+1
)]
This is another implicit method which needs a function solve at each step to nd U
n+1
using a predictor-
corrector method.
It is also a symmetric method i.e. if you know U
n
and you wish to nd U
n+1
with step-size h then
this is the same method for nding U
n
given U
n+1
and step-size h. This property is important for
nding approximations to the solutions equations such as u +u = 0 which are the same both forwards
and backwards in time.
If we now apply this method to u = Au we have
U
n+1
= U
n
+
h
2
[AU
n
+AU
n+1
]
so that
U
n+1
=
_
I
hA
2
_
1
_
I +
hA
2
_
U
n
A straight forward calculation of the eigenvalues of the matrix linking U
n
to U
n+1
shows that growth
rate of the component in the direction of
j
is given by

1 +
h
j
2
1
h
j
2

If we take = p +iq then

1 +
h
2
1
h
2

2
=
_
1 +
ph
2
_
2
+q
2
h
2
4
_
1
ph
2
_
2
+q
2
h
2
4
Thus
< 1 if
_
1 +
ph
2
_
2
+q
2
<
_
1
ph
2
_
2
+q
2
i.e if p < 0
The stability region of the numerical method is thus the half-plane p < 0 illustrated in Figure 6.7
which is identical to the region of stability of the underlying dierential equation. As a consequence the
dynamics of the solution of the trapezium rule exactly mirrors the true dynamics. This is an excellent
state of aairs.
73
The local error LTE of the Trapezium Rule is given by:
LTE =
h
3
|u

12
.
If h is constant the overall error is then proportional to h
2
. To estimate the step-size we keep LTE <
TOL much as before.
The Trapezium Rule, together with a third order method to control the local error, is implemented in
the MATLAB routine ode23t . Here the (second order implicit) Trapezium Rule is a good one to use
if accuracy is not essential. It is very reliable and relatively cheap, but the overall error is quite high
compared to (say) ode45.
30 20 10 0 10 20 30
30
20
10
0
10
20
30
p
q
STABLE
UNSTABLE
Figure 6.7: Stability Region for the Trapezium rule and the Implicit Mid-Point rule.
6.4.5 The Implicit Mid-point Rule...an ideal sti solver?
The implicit mid-point rule is a symmetric Runge-Kutta method closely related to the trapezium rule.
It is given by
U
n+1
= U
n
+hf
_
1
2
(U
n
+U
n+1
)
_
(6.9)
When applied to the linear ODE u

= Au the method in (6.9) gives exactly the same sequence of iterates

as the trapezium rule check this. Thus its stability properties are identical to the Trapezium Rule and
hence are optimal. Like the Trapezium Rule the global error of the Implicit Mid-Point Rule varies as
h
2
and a function solve is required to nd U
n+1
. It has various other very nice features. In particular
it preserves linear and quadratic invariants, so that if u is a solution of the dierential equation and
there exists a vector c and a matrix A so that c u is constant and u
T
Au is a constant then the same
identities hold for the discrete solution as well. It is also a symplectic method (like the St ormer- Verlet
method). An example of the usefulness of these properties comes from the computation of the orbits
of the planets in the solar system. A linear invariant of this system is the linear momentum and a
quadratic invariant is the angular momentum. Both are exactly conserved by this method, which also
comes close to conserving the total energy. Great stu, but at the cost of an expensive function solve.
The Implicit Mid-point rule is the simplest example of a sequence of implicit Runge-Kutta methods
called Gauss-Legendre methods. If you want to solve an ODE very accurately for long times with
excellent stability, but regardless of cost, then these are the methods to use. Gauss-Legendre methods
are the Lamborghinis of the numerical ODE [Link] are expensive and hard to drive, but they have
style! They are implemented in Fortran in the excellent AUTO code.
74
6.4.6 Multi-step methods
These are the most widely used methods for sti problems and include Adams and BDF methods and
the Trapezium rule. There is a vast literature on them, see Iserles. Multi-step methods are very exible
to use and are relatively easy to analyse. They take the form
k

l=0
a
l
U
nl
= h
k

l=0
b
l
f
nl
where f
k
= f (t
k
, U
k
), t
k
= (k 1)h.
To nd U
n
you need to know U
n1
, . . . , U
nk
. To start the method given U
0
you need to use another
method (e.g. RK4) to nd U
1
, . . . U
k1
.
In these methods you use information from previous time steps and (in an explicit method) one function
evaluation to nd U
n
. Thus they are cheaper than RK4 and potentially more accurate. Some properties
of these methods are as follows.
The method has order if the truncation error at each stage is given by Ch
p+1
u
(p+1)
and the overall
method has error of order p, so that the global error is proportional to h
p
.
These methods are implicit if b
o
= 0 and explicit if b
0
= 0. Explicit methods give U
n
directly
whereas implicit methods need an equation to be solved.
They have backwards dierence form if b
o
= 0, b
l
= 0 otherwise
DEFINITION
A multi-step method is A-stable if when used to solve the equation
u

= Au,
when all eigenvalues of A have negative real part, then |U
n
| 0 as n , for all values of h > 0.
More precisely, a method is A-stable if the roots z of the polynomial equation a
l
z
l
hb
l
z
l
= 0
have modulus less than or equal to one if the real part of is less than or equal to zero.
A-stability is a very desirable property for sti problems, which the Trapezium Rule has. However, it
is almost unique in having this property as the following theorem shows:
THEOREM Only implicit methods of order less than or equal to 2 can be A-stable.
Implicit Multi-step methods involve solving non-linear equations. As before the usual method to do
this is a predictor corrector method namely you generate an approximation

U for U
n
using an explicit
predictor method and then correct this to nd U
n
. An easy corrector is to use an iterative one. Suppose
we set U
o
n
=

U
n
where

U
n
is a predicted value obtained using an explicit method. We now perform
the following iteration to nd a sequence of approximations U
r
n
to U
n
:
U
r
n
=
_

l=1
a
l
U
r1
nl
+h

b
l

f
nl
_
/a
o
where

f
k
= f(t
k
, U
r1
k
),
iterating either to convergence or a xed number of times.
75
As in the Runge-Kutta methods, it is common to use two multistep methods simultaneously. One to
perform the numerical solution. The other to give an estimate of the error. In the celebrated Gear
solvers (named after their inventor [Link]) the method chooses what multi-step method to use at each
step from a range of methods of dierent orders (and stability ranges), based on an estimate of the error
from the two methods (the latter is called the Milne device).
MATLAB has an excellent sti Gear solver given by ode15s . It is used in exactly the same way as
ode45 and is the routine to use for most problems, if you dont know much about their structure.
***IF YOU LEARN NOTHING ELSE FROM THIS CHAPTER IT IS TO USE ODE15S
An important sub-set of multi-step methods are BDF (backwards dierence form) methods. BDFn
methods have order n (global errors proportional to h
n
) and BDF1 is just the Backwards Error method.
The important Fortran ODE code DDASSL and ode15s are both based on BDF methods and their
extensions. The rst three BDF methods are given by:
BDF1 : U
n
U
n1
= h f
n
[Backward Euler]
BDF2 : U
n

4
3
U
n1
+
1
3
U
n2
=
2
3
hf
n
BDF3 : U
n

18
11
U
n1
+
9
11
U
n2

2
11
U
n3
=
6
11
h f
n
These methods have excellent stability properties: BDF2 is A-stable (damping out errors but not being
too dissipative). BDF3 is A() stable rather than A-stable; its stability region includes a wedge of angle
and this includes the eigenvalues of many problems such as those arising in uid mechanics.
BDF methods are the Land Rovers of ODE solvers. They may not be pretty or easy to use, but they
will handle rough country and will nearly always get you where you want to go (although they are not
to be used for Hamiltonian problems!).
6.5 Dierential Algebraic Equations (DAEs)
An important application of BDF methods is to dierential algebraic equations. These equations com-
bine algebraic and dierential equations and they are VERY common in many applications for example
chemistry, electronics, uid mechanics and robotics. In a DAE we think of U =
_
p
q
_
as satisfying the
equations
p = f (p, q), 0 = g (p, q). (6.10)
The second of these equations is the algebraic (or constraint) [Link] will rst look at the scalar
case. If we dierentiate this with respect to t we have
0 =
g
p
p +
g
q
q.
If dg/dq is invertible we then have
q =
_
g
q
_
1
g
p
p =
_
g
q
_
1
g
p
f
Thus we arrive at a dierential equation for q. We call problems of this form index-1 problems. If we
have to dierentiate twice with respect to t to get a dierential equation for q we call this an index-2
problem, and if three dierentiations are needed an index-3 problem. Electronics and chemistry tend to
lead to index-1 problems, uid mechanics to index-2 problems and robotics to index-3 problems. MAT-
LAB can handle index-1 problems but has diculties with problems of a higher index. Full details on
76
the theory and computation of DAEs are given in the very readable book by U. Ascher and L. Petzold.
It is not at all obvious how to apply an explicit method to solve such equations, especially to ensure
that the constraint g(p, q) = 0 is met at each stage.
However, it is easy to code these up using a BDF method. For example if we apply BDF2 to the DAE
(6.10) we get:
p
n

4
3
p
n1
+
1
3
p
n2
=
2
3
h f (p
n
, q
n
)
0 = g(p
n
, q
n
)
As before we have to solve Nonlinear equations to nd p
n
and q
n
, but these equations are no worse to
solve than before. In ddassl these equations are solved using quasi-Newton methods in which conjugate-
gradient methods are used to speed up the linear algebra.
The code ode15s can solve DAEs of the form
M u = f(u)
Here the matrix M can be singular, for example if
M =
_
1 0
0 0
_
, f
_
f
g
_
and u
_
p
q
_
then we have p = f and 0 = g as in (6.10).
When using ode15s in this way you specify M in advance and then proceed much as before. More
details are given by help ode15s.
77
Chapter 7
Two Point Boundary Value Problems
(BVPs)
7.1 Introduction
Two point boundary problems (2pt BVPs) take the form:
u = f (x, u) u IR
n
g
1
(u(a)) = x IR
g
2
(u(b)) =
where g
1
, IR
m
and g
2
, IR
nm
. Here x is usually thought of as a spatial variable.
Two point BVPs arise in many contexts of which the following are examples.
1. They are one-dimensional elliptic equations and arise in their own right as descriptions of physical
problems. For example, the equation for the deection of the Euler strut is given by:
d
2
u/dx
2
+sin(u) = 0 u(0) = u(1) = 0,
and the equation for rock deformation (see work by CJB and Professor Giles Hunt) by:
d
4
u/dx
4
+Pd
2
u/dx
2
+f(u) = 0, u(0) = u

(0) = u(1) = u

(1) = 0
2. As steady states of parabolic or hyperbolic partial dierential equations. For example the limit of
the solutions of
u/t = u/x f (x, u)
in the limit of t [Often a good way to solve the steady state problem is to convert it to such
a time-dependent problem].
3. As travelling wave solutions of partial dierential equations or (more generally) as similarity
reductions reductions of partial dierential equations. (These will be described in the semester
Two course on Mathematical Modelling.) An example of this is given by the Bergers equation:

t
+

z
=

z
2
If we look for a travelling wave solution with
(z, t) = u(x) where x = z ct
78
and
(, t) = 1, (, t) = 1.
Then u satises the BVP
c
du
dx
+u
du
dx
=
d
2
u
dx
2
u() = 1, u() = 1.
which is a two-point boundary value problem.
A special (and very hard) case of two-point BVPs is given when a = and/or b = . Here we look
for ground state, homoclinic or heteroclinic solutions. These are very common in PDEs as self-similar
solutions or in quantum mechanics as bound-states We have just seen such a problem in the travelling
wave example.
There is a huge dierence between BVPs and IVPs. In general an initial value problem always has a
unique solution. However a BVP may have one, several (sometimes ) or no solutions.
The theory of such problems is very subtle. An excellent description of this and of numerical methods
for them is given in the book on 2pt BVPs by Ascher, Matheij and Russell. There are many methods
for solving BVPs. These include shooting methods, nite dierence methods, nite element methods,
spectral methods and collocation methods. The latter are used in the MATLAB code bvp4c and
in Fortran codes such as COLSYS, MOVCOL and AUTO. They are especially suitable for use with
adaptive and non-uniform meshes. In this course we will look at shooting methods and nite dierence
methods.
7.2 Shooting Methods
These are the simplest methods for solving BVPs and are based on using IVP software such as ode15s.
They are easy to use and it is often worth trying them rst regardless of the nature of the underlying
problem, although they can perform badly on problems with boundary laters. In a shooting method you
reduce the problem to one of nding the correct boundary values. As this is such a neat idea you should:
SHOOT FIRST AND ASK QUESTIONS LATER!
The idea behind shooting methods is simple. Let y be a solution of the initial value problem
dy
dx
= f (x, y), y(a) = z IR
n
(7.1)
Suppose we take general initial conditions y(a) = z. The initial value problem (7.1) can always be
solved for such general values so that y(b) is then a function F(z) of z. If we can nd a suitable vector
z so that
g
1
(z) = and also g
2
(F(z)) =
then we will have found a solution of the corresponding boundary value problem by simply setting
u(x) y(x).
The procedure works as follows
1. At the point x = a we have g
1
(z) = . As g
1
IR
m
this condition xes m components of the
vector z say z
1
, . . . , z
m
. The other n m components (z
m+1
, . . . , z
n
) can be chosen freely.
2. Next use an initial value solver such as ode15s to calculate y(b) F(z).
79
3. To satisfy the boundary conditions we must have g
2
(F(z)) = . However, as g
2
IR
nm
we have
n m (usually nonlinear) equations to be satised for the n m unknowns z
m+1
, . . . , z
n
.
These nonlinear equations can be solved using a nonlinear solver such as the Newton-Raphson
algorithm or a MATLAB routine such as fzero or fminsearch.
Shooting methods perform well for problems without boundary or internal layers. For those with
boundary layers the eects of exponentially growing terms make these methods hopelessly ill conditioned
and unusable. There are extensions of shooting methods called multi-shooting methods which can be
used for problems with boundary layers. However these are very problem dependent and are not easy
to use.
Example 1
The easiest example of an application of these methods is to linear boundary value problems of the form
a(x)
d
2
u
dx
2
+b(x)
du
dx
+c(x)u = d(x)
with the boundary conditions u(0) = 0, u(1) = 0.
We consider the initial value problem
ay

+by

+cy = d (7.2)
with y(0) = z
1
and y

(0) = z
2
The rst boundary condition on u forces z
1
= 0. However z
2
is arbitrary at this stage.
Because of the linearity of (??) it follows that there are constants p and q so that
y(1) = p +qz
2
[Exercise: prove this]
The values of p and q can be found easily by using an initial value solver applied to (??). For example,
if we set z
2
= 0 then p = y(1) and if z
2
= 1 then q = y(1) p. The value of z
2
for which y(1) = 0 is
then given by
z
2
= p/q.
Example 2
Suppose now that we want to solve the nonlinear boundary value problem
u

+u
2
= 1, u(0) = u(1) = 0.
as before, we solve the initial value problem
y

+y
2
= 1, y(0) = z
1
, y

(0) = z
2
, (7.3)
with z
1
= 0 forced. In this case we have
y(1) = F(z
2
)
where F is a nonlinear function of z
2
. Applying the second boundary condition, y(1) = 0, it follows
that z
2
must satisfy the equation
F(z
2
) = 0.
The function F(z
2
) can be computed by using ode15s. This is illustrated below. It is clear that this
function has a single root at about z
2
= 33.46
80
0 5 10 15 20 25 30 35 40 45 50
8
6
4
2
0
2
4
6
root at 33.4618
z
2

f
The code to generate the function F(z
2
) is given below:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Shooting code %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=funn(x)
[t,u] = ode15s(u2,[0 1],[0 x]);
s = size(t);
ss = s(1);
F = u(ss,1);
--------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% u + u^2 = 1 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f=u2(t,u)
f=zeros(2,1);
f(1) = u(2);
f(2) = 1-u(1)^2;
--------------------------------------------------------
One way to nd z
2
is to use the MATLAB command fzero given an initial guess of z
2
= 30. Notice
the way that this code uses the function funn computerd using ode15s.
>> z2 = fzero(funn,30)
z2 =
81
33.4618
The resulting solution u and u

is illustrated below:
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
40
30
20
10
0
10
20
30
40
x
u
du/dx
An alternative method, which exploits the structure of the problem, is to use a Newton-Raphson method.
As F(z
2
) = y(1) it follows that
dF
dz
2
=
dy(1)
dz
2
= (z
2
)
where the function (x) satises the variational equation

+ 2y = 0, (0) = 0,

(0) = 1.
To use this method we take a guess z
(n)
2
for z
2
and then simultaneously solve the two ODES
y

+y
2
= 1, y(0) = 0, y

(0) = z
(n)
2
,

+ 2y = 0, (0) = 0,

(0) = 1,
using ode15s. We then use a Newton-Raphson iteration to improve the guess via the equation
z
(n+1)
2
= z
(n)
2

y(1)
(1)
.
This procedure converges rapidly given a reasonable starting guess.
7.3 Finite Dierence Methods
7.3.1 Discretising the second derivative.
Finite dierence methods aim to to nd the solution at all interior points within the interval [a, b].
Suppose that we divide up the interval [a, b] into N equal divisions, so that if h = (b a)/(N 1) we
set x
i
= a + (i 1)h i = 1, . . . , N.
Now we introduce a vector U
i
so that
U
i
u(x
i
) (7.4)
Now, by using a Taylor series expansion, it follows that
u(x
i+1
) = u(x
i
) +hu

(x
i
) +
h
2
2
u

(x
i
) +. . .
u(x
i1
) = u(x
i
) hu

(x
i
) +
h
2
2
u

(x
i
) +. . .
By making a slightly more detailed calculation it follows that
u(x
i+1
) 2u(x
i
) +u(x
i1
)
h
2
= u

(x
i
) +
h
2
12
u
iv
(x
i
) +. . . (7.5)
82
We now can use (??) to approximate u

by
u

U
i+1
2U
i
+U
i1
h
2


2
h
2
U
i
(7.6)
where
2
is the iterated central dierence operator. The leading error made in (??) when approximating
u

by the iterated central dierence is then proportional to h

2
u
iv
. Now suppose that A is the tri-diagonal
matrix given by
A =
_

_
1 0 . . . 0
1/h
2
2/h
2
1/h
2
.
.
.
.
.
.
.
.
.
1/h
2
2/h
2
1/h
2
0 1
_

_
,
If the vector U then satises the linear system
AU =
_

f
2
.
.
.
f
N1

_
where f
i
f(x
i
) (7.7)
then (??) is a discretisation to the second order BVP with Dirichlet boundary conditions which is given
by
u

= f(x), u(a) = , u(b) = .

If we replace U by u in (??) the expression (??) implies that we would make an error of order h
2
as
h 0. Note that the structure of A forces U
1
= and U
N
= . The vector U can then be found via
the operation
U = A
1
_

f
2
f
N1

_
Provided that the problem is not too irregular it can then be shown that in the limit of N (h 0)
we have
U
i
= u(x
i
) +O(h
2
) as h 0.
The solution of the BVP thus becomes (under discretisation) the problem of solving the linear system
(??).
If f = [, f
2
, . . . , f
N1
, ]

then this can be done in MATLAB via the command A\f.

Whilst this is simple to use it is not especially ecient as the backslash operator does not exploit the
special tri-diagonal structure of Matrix A. A much more ecient algorithm of complexity O(N) is the
Thomas algorithm which exploits this structure and is based on the LU decomposition. This algorithm
is described in Assignment 4. (Note, use of the MATLAB sparse matrix routines will speed things up
here.)
Eciency is important as the relatively large error of O(N
2
) in this method means that we have to
take a large value of N to get a reasonable error estimate. More careful discretisations lead to smaller
errors at the expense of more work.
7.3.2 Dierent boundary conditions
We can extend this idea to BVPs with Neumann boundary conditions. For example the BVP u

= f
with the boundary conditions
u

(a) = , u(b) = .
83
Problems such as this arise frequently in models of heat conduction and electrical ow. For example,
the boundary condition u

(a) = 0 corresponds to a thermal or electrical insulator. There are many ways

to deal with such a derivative boundary condition. In the simplest we approximate
u

(a) by (U
2
U
1
)/h (7.8)
This will change the matrix A to the tridiagonal matrix
A
_

_
1/h 1/h 0 . . . 0
1/h
2
2/h
2
1/h
2
.
.
.
.
.
.
.
.
.
1/h
2
2/h
2
1/h
2
0 0 1
_

_
and we must consider the solution to the linear system:
AU =
_

f
2
.
.
.
f
N1

_
(7.9)
The linear equation (??) is then a discretisation of the BVP
u

= f u

(a) = , u(b) = .
The expression (??) is only accurate to O(h) as an approximation to u

(a). We can improve this by

using the approximation
u

(a) (4U
2
U
3
3U
1
)/2h
Exercise. Show that this expression is accurate to O(h
2
)
This leads to a slightly dierent matrix A. Unfortunately the resulting matrix is no longer tri-diagonal
and we cannot use the Thomas algorithm to invert it.
A special case of Neumann boundary conditions arises when u

(a) = 0. In this case we can introduce

a ghost point U
0
approximating u(a h). The boundary condition implies that to 0(h
3
) we have
U
0
= U
2
. At the point x = 0 we have
u

(a)
U
0
+U
2
2U
1
h
2
=
2U
2
2U
1
h
2
(7.10)
This approximation to u

(a) is correct to O(h

2
). (Can you show this?)
The resulting matrix A is then given by
A =
1
h
2
_

_
2 2
1 2 1
1 2 1
.
.
.
.
.
.
_

_
Note that in this case A is tri-diagonal, and we can again solve the linear system by using the Thomas
algorithm.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84
% %
% Code to solve the Neumann two-point boundary value %
% problem: %
% %
% u = -exp(x), u(0) = 0, u(1) = 0. %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set up the mesh
%
N = 101;
h = 1/(N-1);
x = [0:h:1];
%
% Set up the matrix
%
A = diag(ones(N-1,1),-1) + diag(ones(N-1,1),1) - 2*diag(ones(N,1));
A = A/h^2;
%
% Neumann boundary condition at x = 0
%
A(1,2) = 2/h^2;
%
% Dirichlet boundary condition at x = 1
%
A(N,N-1) = 0;
A(N,N) = 1;
%
% Set up the right hand side
%
f = -exp(x);
%
% Modify this to allow for the Dirichlet boundary condition
%
f(N) = 0;
85
%
% Solve the system (without using the Thomas algorithm)
%
U = A\f;
%
% Plot the solution
%
plot(x,U)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
x
U(x)
Neumann condition
Dirichlet condition
A second example of dierential equations occurs for BVPs with periodic boundary conditions for which
u(a) = u(b) and u

(a) = u

(b). Periodic boundary conditions often arise when looking at those BVPs
that describe the travelling wave solutions of hyperbolic equations. They also arise naturally when
solving BVPs on circles or spheres. A natural example of the latter being weather forecasting on the
whole globe. Often BVPs with periodic boundary conditions are solved using spectral methods in which
the solution is expressed as a combination of trigonometric functions. However they can also be solved
by using nite dierence methods. Exploiting periodicity we have
u

(a)
u(a +h) +u(a h) 2u(a)
h
2
=
u(a +h) +u(b h) 2u(a)
h
2
(7.11)
which we approximate by
u

(a)
U
2
+U
N1
2U
1
h
2
Setting U = [U
1
, . . . , U
N1
]

(noting that U
1
= U
N
) the resulting matrix A is given by
A =
_

_
2/h
2
1/h
2
0 . . . 0 1/h
2
1/h
2
.
.
.
.
.
. 0
0
.
.
.
.
.
.
1/h
2
. . . . . . 1/h
2
2/h
2
_

_
The resulting discretisation of u

then has an error of O(h

2
) as before. The matrix A is not tri-diagonal,
but its special periodic structure means that it can be inverted quickly by using the FFT.
86
7.3.3 Adding in convective terms
In many applications we meet convection- diusion problems which typically take the form
u

+u

= f with small.
These arise in uid mechanics problems, where is a measure of the viscosity of the uid. The most
accurate discretisation of u

is given by the central dierence discretisation

u

U
2h

U
i+1
U
i1
2h
A simple calculation shows that the central dierence gives an error O(h
2
) when approximating u

.
However this discretisation of u

leads to problems. In particular A

1
may be nearly singular and
ill-conditioned. If is zero then A has an eigenvector of the form e = (1, 1, 1, 1, 1, 1, . . .) with zero
eigenvalue eigenvalue. Check thisso that A
1
is singular.
For small , A continues to have a similar eigenvector to with a small eigenvalue. When inverting A
the contribution of e to the solution is greatly amplied, and can lead to an oscillatory contribution to
the solution. This instability is BAD but can easily be detected as it is on a grid scale. i.e all oscillations
in the solution are the same size as the mesh.
Oscillations can be avoided by taking h small enough to avoid this instability so that
h 2
However, this is a very restrictive condition if is small.
One way to avoid the instability when > 0 is to replace u

by its upwind dierence

U
i+1
U
i
h
This is a less accurate approximation of u

than the central dierenceas

U
i+1
U
i
h
2
= u

+O(h)
However the resulting matrix A is much better conditioned and inverting it does not lead to spurious
oscillations in the solution What this approximation is doing is to exploit a natural ow of information
in the system.
This is not always a good solution as often we cant tell in advance if we need to use the upwind dier-
ence or the downwind dierence
U
i
U
i1
h
. However, like a sti solver, the use of the upwind dierence
allows us to use a step size governed by accuracy rather than stability. An example of the eect of using
these dierent discretisations is given in Assignment 4.
7.4 Nonlinear Problems.
So far we have looked at linear two point BVPs. However,
MOST PROBLEMS ARE NONLINEAR
87
An important class of nonlinear two point BVPs (called semilinear problems) take the form
a(x)u

+b(x)u

+f(x, u) = 0, (7.12)
and we need to develop techniques to solve them. Two examples of such problems are given by the
dierential equation
u

+e
u/(1+u)
= 0, u(a) = u(b) = 0
which models combustion in a chemically reacting material, and
u

+
2
x
u

+k(x)u
5
= 0, u

(0) = 0, u() = 0, u > 0,

which describes the curvature of space by a spherical body.
A special case is given by the equation
u

+f(u) = 0, u(a) = , u(b) = . (7.13)

Discretising the dierential equation (??) leads to a set of nonlinear equations of the form
AU +f(U) = 0 (7.14)
where A is the matrix given earlier and the vector f is given by
f
i
f(U
i
).
As with most nonlinear problems, we solve this system by iteration starting with an initial guess
although, be warned, the system may have one, none or many solutions. Perhaps the most eective
such method is the Newton-Raphson algorithm. Suppose that U
(n)
is an approximate solution of (??),
we dene the residual R
(n)
by
AU
(n)
+f(U
(n)
) = R
(n)
. (7.15)
The size of R
n
is a measure of the quality of this solution.
Now if we dene the Jacobian of the nonlinear function f by
J
ij
= f
i
/U
j
the linearisation L of (??) acting on a vector is given by
L A +J.
The Newton-Raphson iteration updates U
(n)
to U
(n+1)
via the iteration
U
(n+1)
= U
(n)
L
1
R
(n)
Here the vector W L
1
R
(n)
can be found by solving the linear system
AW +JW = R
(n)
.
This algorithm converges rapidly (often in about 5 iterations) if the initial guess U
(0)
is close to the true
solution. Finding such an initial guess can be dicult for a general problem and often requires some
a-priori knowledge of the solution.
The Fortran code AUTO uses the Newton-Raphson method coupled to a path-following (homotopy)
method to nd the solution of a version of the nonlinear systems arising from discretisations of nonlinear
BVPs obtained by using collocation.
88
7.5 Other methods
At present collocation is the best method for solving 2pt BVPs. This is a bit like the Finite dierence
method but uses much higher order polynominals. It is also closely related to implicit Runge-Kutta
methods. However, the collation method has no easy extension to higher dimensions. In higher dimen-
sions two-point BVPs generalise to elliptic partial dierential equations. For such problems the nite
element method is the mainly used solution procedure, although spectral methods are widely used for
problems with a relatively simple geometry e.g. weather forecasting on the sphere.
CJB 2006
89