Math 563 Lecture Notes

Numerical methods for boundary value problems

Jeffrey Wong
April 12, 2020

Related reading: Ascher and Petzold, Chapter 6 (a good discussion of stability ) and
Chapter 7 (which includes details on multiple shooting and setting up Newton’s method for
these problems).

1 Boundary value problems (background)

An ODE boundary value problem consists of an ODE in some interval [a, b] and a set of
‘boundary conditions’ involving the data at both endpoints.

After converting to a first order system, any BVP can be written as a system of m-equations
for a solution y(x) : R → Rm satisfying
dy
= F (x, y), x ∈ [a, b]
dx
with boundary conditions
B(y(a), y(b)) = →0.
As with ODE IVPs, this ‘standard’ form is typically used for numerical computation (e.g.
Matlab’s bvp4c requires the ODE written this way), although there are often simplifications
that can be made given some structure.

For example, the two-point boundary value problem (i.e. a BVP with a second order
ODE and two boundary conditions)

yrr = 4y − yr, x ∈ [a, b], y(a) = c, yr(b) = d

can be written in this form, where w = (y1, y2) with w1 = y and w2 = yr and

r
= w2
w1r w1(a) −
w2 = 4w1 − w2 B(y(a), y(b)) =
c w2(b) −
d

1
General boundary conditions can be difficult to deal with. The most common (and simplest)
type are linear boundary conditions, where

Bay(a) + Bby(b) = c, Ba, Bb = m × m matrices.

The boundary conditions are separated if they can be written as

Bay(a) = ca, Bby(b) = cb

e.g. for the example above (where Ba = [ 1 0 ] and Bb = [ 0 0 ]).
00 01

For an IVP, the effect of the initial condition is carried forward by the ODE from the
initial point. We can follow this data and take steps one at a time, which is natural and
straightforward.

In contrast, a BVP is much harder to solve because the ODE carries information from
the boundaries throughout the interval. We cannot start at one point and step forward,
because that would ignore the other boundary data. Thus, more work is required.

1.1 Stability
The stability of BVPs is a substantial topic; a quick study will help to illustrate the issue.
First, suppose we have a linear IVP

yr = Ay + p(t), y(0) = y0

where A is an m × m matrix. If the ODE function is perturbed slightly,

y˜r = Ay˜ + p˜(t), y(0) =

0 then the difference u = y˜ − y evolves according to

ur = Au + δ(t), u(0) = 0, δ := p˜(t) − p(t).

From ODE theory, we know that the eigenvalues of A determine stability in the usual way
(stable if Re(λ) < 0). A similar result holds for changes in the initial condition, and a linear
stability analysis (as we did before) extends this to the non-linear case.

Now contrast with the BVP case. Consider a linear, separated BVP

yr = Ay + p(t), Bay(a) = ca, Bby(b) = cb

and a perturbed problem (assuming the BCs do not change for simplicity)

y˜ r = Ay˜ + p˜(t), Ba y˜(a) = ca, Bb y˜(b) = cb.

Then the difference u = y˜ − y solves

ur = Au + δ(x), Bau(a) = Bbu(b) = 0.

2
That is, the (linear) stability of the BVP depends on the properties of a boundary value
problem. The solution can be written in terms of a Green’s function:
∫ b
u(x) = a (p˜(y) − p(y))G(x, y)
dy
and thus stability really depends on the nature of the Green’s function (which we will not
discuss here!). In particular,
∫ b
u ≤ G |p˜ − p| dy, G ∞ := max |G(x, y)|
a a≤x,y≤b
∞

which describes the sensitivity of the solution to changes in the ODE function.1

The point here is that the notion of stability from IVPs, e.g. that yr = λy is stable as
t increases if Re(λ) < 0, does not qutie apply here. The boundary conditions can both
create new instability or prevent it.

As a trivial but useful example, consider, for c > 0,

1 2
r −c0 0 c
y = Ay, A= , y (0) = 1, y (b) = 1

This is nothing more than two uncoupled IVPs:

yr1 = −cy1, y1(0) = 1

r
y2 = cy2, y2(b) = 1.
The equation for y1 is solved forward (from t = 0 up to t = 1), and is stable since Re(c) > 0;
the result for IVPs apply here since it does not care about the boundary condition on the
other side.

The equation for y2 is solved backward (down to t = 0); it is equivalent to the IVP
dy2 = −cy , y (τ = 0) = 1
2 2
dτ
for τ = b − t, and so it is also stable. It follows that the BVP must be stable as well (by
any reasonable definition).

However, an initial value problem for the system,

yr = Ay, y(0) = y0
1
Note: If the BCs are also allowed to be perturbed, then we instead get

u ≤ (max{Φ , G })(change in BCs +∞ where Φ is a fundamental

3
 ∫ b
matrix for y′ = Ay. |p˜ − p| dy)
a

4
is clearly not stable since one component grows exponentially (and this is also true in the
other direction). The ‘instability’ of the IVP does not apply to the BVP because the bound-
ary stops any unstable growth.

The point: The key practical result here is that the stability of the boundary value problem
and the associated initial value problem (in whatever direction) are not necessarily the
same. Thus, using ‘IVP’ techniques to solve the boundary value problems can introduce
new (artificial) instability [we will see this shortly].

The second point is that the ‘stiffness’ condition for a BVP is different. The linear,
constant coefficient BVP
yr = λy, x ∈ [0, b] ...BCs...
may be considered stiff if
b|Re(λ)| 1.
The condition here generalizes as in the IVP case. This means that a solution y = eλx
may grow rapidly relative to the interval size (plug in x = b). Because there is no longer
a ‘forward’ direction as in an IVP, stiffness does not require Re(λ) < 0.

Stiff BVPs can be challenging to solve because they often involve boundary layers, thin
regions of rapid variation near the boundaries (or elsewhere).

1.2 Uniqueness
One other subtlety is that existence/uniqueness is more complicated for BVPs.

Recall that for an IVP

y = F (t, y), y(0) = y0
there is a unique (local) solution under mild assumptions - namely, that F is continuous and
Lipschitz in y. In contrast, a BVP may have:
i) A unique solution
ii) multiple solutions
iii no solution at all
Obviously, a problem with no solution cannot be solved numerically. Case (iii) also has
consequences for solvable problems - it means that even if a problem has a solution,
nearby problems (which we may encounter in the approximation process) might not.

From a numerical perspective, (ii) is not a serious problem, since numerical methods will
naturally select one solution in the computation. We just need to be careful to find the right
solution - for instance, by using an initial guess close to the desired solution.

One important example is an eigenvalue problem such as

yrr = −λy, y(0) = y(π) = 0

5
which has solutions y(x) = sin(kx) for λ = k an integer and no solutions otherwise. Note
that the eigenvalue problem can be regarded as a system of ODEs; set w1 = y and w2 = yr
and add a third trivial variable λ:
r
w 1 = w 2, wr2 = −λw1, λr = 0.

This is a first order system of three ODEs for (w1, w2, λ). The trick of replacing a param-
eter with a trivial ODE is a ‘reformulation trick’, which we can use to solve for unknown
constants.2

2 Boundary value problems (shooting, part I)

To start, we consider a typical two-point boundary value problem

yrr = f (x, y, yr), x ∈ [a, b], y(a) = c, y(b) = d

for a function y(x). Unlike an initial value problem, there are conditions involving y at both
endpoints of the interval, so we cannot just start at x = a and integrate up to x = b.

One approach (to be considered later) is to approximate y(x) by some finite difference or
other approximation scheme, then arrive at a system for the discrete values y(xi) that can
be solved all at once. The downside is that one needs to typically solve large linear systems
that require some effort.

However, there is a cheaper way. Note that y(a) is specified. We can use an IVP solver
to compute a solution y(x; s) given initial conditions

y(a; s) = c, yr(a; s) = s

where s is a ‘guess’ at the correct value of yr(a). This solution will hit x = b at some value
that is probably not equal to the correct value, y(b) = d. However, we can then adjust the
value of s, refining the guess until is is (nearly equal to) the right value.

This method is called shooting - the analogy here is to shooting a ball [insert relevant
sport here] at a goal, determining the unknown correct velocity by throwing it too fast/too
slow until it hits the goal exactly.

In particular, it means that if you have a routine to solve a BVP in standard form, you can also solve
2

BVPs with unknown parameters with the same code by adding trivial equations.

6
 2.1 The details
Here are the details for using shooting to solve the two-point BVP

yrr = f (x, y, yr), x ∈ [a, b], y(a) = c, y(b) = d. (BVP)

1) Setup the IVP: First, we set up the shooting initial value problem

yrr = f (x, y, yr), y(a) = c, yr(a) = s. (IVPs)

The solution to (BVP) is a solution to (IVPs) for a certain value s∗ of s that we must find.
Let us define
y(x; s) = solution to (IVPs) for that value of s.
2) Define the goal: Next, we define the ‘goal function’ that tells us what the correct value
s∗ must satisfy in the form G(s) = 0. Here, the goal function is

g(s) = y(x; s) − d.
x=b

We seek s∗ such that g(s∗) = 0. Thus, the problem has been reduced to a non-linear
equation that can be solved via the method of our choice. Note that computing g(s)
involves solving (IVPs), so each evaluation of g is fairly expensive. Newton’s method is then
a desirable method due to its fast convergence.

2b) Setup variational problem for Newton: If using a ‘derivative free’ method like
the secant method, this step can be skipped. 3 To use Newton’s method, we also need the
derivative of g. This requires knowing the derivative of y with respect to s. Let
∂y(x; s)
z(x; s) = .
∂s
Then, by differentiating the DE and initial conditions, we get the variational ODE
∂f ∂f r
zrr = z+ z , z(a; s) = 0, zr(a; s) = 1 (IVPZ )
s
∂y ∂yr
which describes the change in the IVP solution with respect to its parameter s.

The derivative of g involves z at the endpoint x = b:

gr(s) = z(b; s).

Thus, computing gr(s) requires solving the IVP (IVPZs). Putting this together, we now can
compute g(s) and gr(s) so one just has to use Newton’s method:

s =s
g(sn)
n+1 —
n
gr(sn)
3
There are two reasons to do so: First, if the variational ODE is complicated or not smooth enough and the
advantages of Newton’s method is not needed, a method that needs to know only the ODE function is desir-
able. Second, the variational ODE might be ill-behaved, introducing more possibility of error/instability/etc.

7
which will hopefully converge to the correct value s∗ if the initial value is good enough.
Explicitly, the iteration reads

s =s
y(b; sn) − d
n+1 n — ∂y .
∂s
(b; sn)

To summarize, when using Newton’s method, the steps are:

• Setup the shooting problem and define the goal function g(s)

• Use Newton’s method to solve g(s) = 0 for the ‘correct’ value s∗. At each step,
-solve the shooting IVP (IVPs) to compute g
-solve the variational ODE (IVPZs) on the same grid and obtain gr
Take a step of Newton using g and gr

• Done; the BVP solution is the solution to the IVP with the computed value of s∗

Note that the variational ODE for z(x; s) is solved on the same grid because it must know
the IVP solution y(x; s) at the grid points. You could also just put the two systems together
and solve them at once.

Also note that shooting really only requires

• an initial value problem solver (e.g. RKF)

• a non linear system solver (e.g. Newton)

and any method of each type can be used. The appropriate choice will vary by problem.

2.2 Two-point BVP example

The two-point boundary value problem

yrr = 2y3 − 6y − 2x3, y(1) = 2, y(2) = 5/2 (2.1)

has a solution y(x) that happens to be given exactly by

1
y(x) = x + , x ∈ [1, 2].
x
Let y(x; s) be the solution to the shooting IVP

yrr = 2y3 − 6y − 2x3, y(1) = 2, yr(1) = s. (2.2)

The goal function is

g(s) = y(2; s) − 5/2 (2.3)

8
so that the solution to the IVP with s such that g(s) = 0 is also a solution to the BVP (2.1).

For computation, set w1 = y, w2 = yr and z1 = ∂y/∂s and z2 = ∂yr/∂s. Then

w1r = w2
w2r = 2w31 − 6w1 − 2x3
z1r = z2
z2r = (6w21 − 6)z1
with initial condition
(w1, w2, z1, z2) = (0, s, 0, 1).
The goal function and its derivative are then
g(s) = w1(2; s) − 5/2, gr(s) = w3(2; s).
Note that the first two equations can be solved first, and the variational ODE (the last two
equations) can be solved after if the same grid is used.

The ODE has the feature that if yr(1) is too large (e.g. yr(1) > 1/2) then the solution
to the IVP (2.2) has a singularity in [1, 2]. Thus, it is necessary to choose a starting
value s < 1/2, and must be close enough to 0 that it never enters this singular region.
Or, some other, safer method could be used (e.g. bisection).

2.3 Perils
It is important to note that the stability of shooting depends on the stability of the
IVP, even though (as discussed) the actual stability of the BVP may be different.

In particular, a stable BVP may have a shooting IVP that is very stiff or otherwise ill-
behaved. That is, by using shooting we may introduce numerical issues hat are not
inherent to the problem, which is a significant disadvantage.

Moreover, note that the IVPs must be solved for guesses s that are not correct, so we
may need to solve nearby problems that are worse than the true solution - e.g. if the IVP is
singular for certain values of s.

There are two cheap ways to ‘fix’ shooting if the IVP is not well behaved.

First, if the IVP blows up only in one direction, we can try to shoot from the other side.
[The application is left as an exercise].

A stronger technique is multiple shooting. The idea is that if the IVP wants to blow
up after some distance < b for a BVP in the interval [0, b], we can ‘restart’ the solve by
guessing values inside the interval and then connect the pieces. For the two-point problem
yrr = f (x, y, yr), y(a) = c, y(b) = d

9
multiple shooting with M segments proceeds by picking points x1, x2, · · · , xM—1 in (0, b)
(with x0 = a) and providing a guess at each point. We solve the IVPs

yrr = f (x, y, yr), yr(xj) = sj

for a solution y(x; c), dependent on the guess →s = (s0, · · · sM—1) for the starting value for yr
at each point. The condition on y is a continuity condition; the values y(xj; s) are of course
available if we solve the IVPs in order (j = 0, 1, · · · .)

The goal function is the same as before, g(s) = y(b; s) − d, but now depends on an M -
dimensional vector. The non-linear system can be solved via the usual techniques (e.g.
Newton’s method), but takes some work.

For this reason, multiple shooting is most useful when shooting almost works and it is
reasonably easy to set it up (requiring a small M , or some nice structure). Otherwise, other
techniques (e.g. finite difference methods, to be discussed) should be used. Note that if the
IVP is extremely poorly behaved, then M has to be quite large, which is undesirable as the
size of the system g(s) = 0 to be solved puts a serious demand on the root-finder.

2.4 Continuation
A boundary value problem encountered in plasma physics4 is Troesch’s problem

yrr = µ sinh(µy), y(0) = 0, y(1) = 1

where µ > 0 is a parameter. The IVP for shooting and variational ODE for z = ∂y/∂s are

yrr = µ sinh(µy), y(0) = 0, yr(0) = s

zrr = −µ2 cosh(µy)z, z(0) = 0, zr(0) = 1
and the goal function is
g(s) = y(1; s) − 1.
It can be shown that

yr(0) > 8e—µ =⇒ the IVP solution has a singularity in [0, 1].

Note that this implies that the correct value of s decays at least like e—µ as µ increases,
so the right value can become quite small (this means care is required in using relative vs.
absolute error tolerances in Newton’s method).

More importantly, when µ is not small (about µ > 5), the ODE is sensitive to the value of
yr(0), and one has to pick a very good initial guess to get Newton’s method to converge.

If shooting is used, this problem is a good candidate for continuation; the problem is easy
to solve for µ = 0, and one can increase µ from there.
4
See https://doi.org/10.1016/0021-9991(73)90165-4.

10
3 Non-scalar problems
Consider a boundary value problem

yr = F (x, y), x ∈ [a, b], B(y(a), y(b)) = 0

with a first order system of m equations. The shooting IVP is

yr = F (x, y), x ∈ [a, b], y(a) = s

which has a solution y(x; s). Note that the guess is, at worst, a vector in Rm. If some
components at x = a are known, we can reduce the dimension (e.g. if y1(a) = 0).

The goal function is then

G(s) = B(y(a; s), y(b; s)), G : Rm → Rm

again noting that its dimension may be less if it can be simplified.

We could use a derivative-free method to solve G(s) = 0 (and be done); the details of
such methods will not be discussed here. To use Newton’s method, we need the Jacobian of
G. Define the matrix-valued function (the variation of the solution with respect to s)

Z(x; s) = ∂y(x; s)
.
∂s
That is, Z has columns [z1, · · · , zm] where

zj = ∂y(x; s)
= variation of y with respect to sj.
∂sj
By differentiating the IVP, we see that Z solves the variational ODE

Zr = J(x, y) Z, Z(a; s) = I (3.1)

where J(x, y) = ∂F/∂y is the Jacobian of F with respect to y. Note that the ODE (3.1) is
really a collection of m equations for

zj(x; s) = ∂y(a; s)
, j = 1, · · · , m
∂sj
which is how they would be solve numerically. Because the ODE is linear and the Jacobian
is shared by all m equations, all the zj’s can be computed efficiently at once.

With Z known, we can then calculate the Jacobian of the goal function needed for shooting.

11
3.1 Period finding (VDP)
The Van der Pol equation is (in system form)
r
x1 = x2, x2 r = µ(1 − x12)x2 − x1

which describes a non-linear oscillator with negative damping when its amplitude is small.

For certain values of µ (when µ > 0) this system has a limit cycle - a unique periodic
solution. This can be formulated as a boundary value problem. Let the period be T , and
suppose that x1(0) = 0 (it is not hard to show that the limit cycle must cross the x2-axis,
so we can pick that crossing point as the t = 0 point as the system is autonomous).

The BVP to solve is

x1(0) = 0
xr1r = x2 x1(0) = x1(T )
x = µ(1 − x2)x − x
2 1 2 x (0) = x (T )
1
2 2

where T is also an unknown. We can convert to a nicer form using two tricks. First, scale to
get a fixed domain; set τ = t/T as a new time variable. Second, regard T as an independent
variable with T r = 0 to get a first order system for (x1, x2, T ). The result is

dx1
dτ
dx2 = Tx2
x1(0) = 0
= T (µ(1 − x2)x — x ) x (0) = x (1)
1 2 1 1 1
dτ
dT x2(0) = x2(1)
=0
dτ
which is a first-order system in a standard form (with three BCs as expected). The fact that
T r = 0 means that the process can be simplified a bit.

To apply shooting to solve this BVP, we solve an IVP of the form

dx
= TF (x), x(0) = (0, s) with a guess (s, T )
dτ
noting that the T r = 0 equation doesn’t actually need to be included in the ODE as there is
nothing to solve. The goal function is
x1(1; s)
G(s, T ) = .
x2(1; s) − s
If Newton’as method is to be used, we need the appropriate variational ODEs. Let

12
 ∂x1 ∂x2
z= , w= .
∂s ∂T

13
Then, with J(x) the Jacobian of F (x),
dz dw
= TJ(x)z, = F (x) + TJ(x)w
dτ dτ
with initial conditions

z1(0) = 0, z2(0) = 1, w1(0) = w2(0) = 0.

The Jacobian of G can then be computed using z and w at t = 1, e.g.

∂G2 ∂
= (x (1; s) − s) = z (1; s).
2 2
∂s ∂s
Note that one could instead think of the system as a three-equation system xr = H(x)
for (x1, x2, x3) with xr3 = 0 (as T ) and with a shooting parameter s = (s1, s2) and initial
conditions
x1(0) = 0, x2(0) = s1, x3(0) = s2.
Then the variational ODE for the matrix Z = ∂x/∂s is
∂H
Zr = Z.
∂x
The result is the same as the calculation above after simplifying (we have essentially solved
the xr3 = 0 equation and split Z into its columns).

The shooting method has trouble when the ODE starts to become harder to solve; as µ
increases, the initial guess must be very accurate (as well as the ODE solution). Continua-
tion can help, but one may be better off using another method entirely.

14
4 Finite differences
Suppose, for example, we wish to solve the boundary value problem

yrr = f (x, y, yr), x ∈ [0, b]

(BCs at x = 0 and x = b) (4.1)

This problem can be solved by the following direct process:

1) Discretize the domain: Choose grid points {xn} for the approximation {un}

2) Discretize the ODE: Approximate derivatives with finite differences in the interior
of the domain (where BCs do not apply).

3) Boundary conditions: Discretize the boundary conditions, and modify (2) as needed.

4) Solve the system of equations from (2) and (3) by the usual methods.

Unlike shooting, the entire solution is solved for at once, which avoids the potential insta-
bilities in the auxiliary IVPs. The details depend, of course, on the choice of grid points
(uniform or non-uniform) and the discretization (which formulas to use).

4.1 The basic approach

We solve the example problem (4.1) with ‘Dirichlet’ BCs

y(0) = c, y(b) = d (4.2)

using a uniform grid

xn = hn, n = 0, 1, · · · , N, h := b/N.
Since the values of y are known at n = 0 and n = N , our unknowns are (y1, · · · , yN—1).

Using the standard centered difference formula for yrr and yr and replacing y(xn) with un
(the approximation to be found) gives
un+1 − 2un + un— un+1 − un—1
=f , un , ), n = 1, · · · , N − 1, (4.3)
1
(xn 2h
h2
and applying the boundary conditions gives (trivially)

u0 = c, uN = d.

This is a system of N + 1 equations for the N + 1 values of {un} on the grid, or alternatively
a system of N − 1 equations for the ‘interior’ values (u1, · · · , uN—1). The latter is preferable
for computation. The result is a system of equations for u = (u1, · · · , uN—1),

G(u) = 0, G : R N —1 → R N —1 .

15
The equations are obtained by taking (4.5) and ‘plugging in’ for the boundary values u0, uN :

G (u) = u − 2u u2 − c
+ c − h2f (x , u , ),
1 2
1 1
1 2h
un+1 − un—1
Gn(u) = , ), n = 2, · · · , N −
n+1 — 2un + un—1 — h f (xn , un
2
u 2, 2h

GN —
(u) = d − 2u N — + uN — — h2f (x N — , uN — , d − uN—2 ),
1
1 2
2h 1 1

This system can then be solved using Newton’s method. Note that the Jacobian J = ∂G/∂u
is tridiagonal, i.e. has non-zero entries only on the center diagonal and one above/below;
the linear system for Newton steps can then be solved efficiently (each linear solve
requires O(N ) operations).

Note that if the BVP is linear then the system has the form Au = →b, so it can be solved
as a linear system. For example, if f = x2 (so yrr = x2) then

 −2 1 0 ···  −c  2 
 x1
0 0  
1 −2 1 ···
 0
    x22 
 . . . 
 2
..
.
. , →
b=
. +h .. 
. . .
A=
    2 

0 0 · .· · 1. −2 1  0 x N — 
0 ··· 0 1 −2 −d x2N —
2
1

Note that the boundary conditions add to the right hand side - every term that does not
depend on u gets moved to the RHS (b) and every term involving u is on the LHS (Au).

4.2 Example (linear)

The Airy function Ai(x) is a multiple of5 the unique solution to

yrr = xy, y(0) = 1, lim y(x) = 0

x→∞

Because the ODE has another solution Bi(x) that grows rapidly (see homework 6), using
shooting can be difficult (it is sensitive to the choice of yr(0), and solutions do not want to
go to zero).

We solve the BVP instead, replacing infinity with a large value L. Letting

xn = nh, h = L/N,

we discretize as discussed above to get

un+1 − 2un + un—1 = h2xnun, n = 1, · · · , N − 1

16
5
to be precise, the Airy function has Ai(0) = 1/(32/3Γ(2/3)).

17
Then the system to solve (which is linear) is Au = b where
 
 
−2 − h2x1 1 0 ··· 0 −1
1 −2 − h2x2 1 ··· 0 0

A= 0
..
.
..
.
..
.  →b = . 
. ,
.
.
   
 0 ··· 1 −2 − h2xN—2 1  0
0 ··· 0 1 −2 − h2xN—1 0

The method has no trouble computing a solution. Note that because x > 0, the matrix A is
diagonally dominant (the diagonal entry is larger in size than the absolute sum of other
entries in the row). It follows that the matrix is invertible and that the system can be stably
solved using LU decomposition without pivoting.

4.3 Neumann (no known values)

If a derivative is specified (‘Neumann boundary conditions’), such as

yr(0) = c, yr(b) = d (4.4)

there is one subtlety. We set up the grid as before and discretize the ODE in the same way
for interior points.
However, no values of y are known, so all N + 1 values of {un} (from n = 0 to n = N )
are unknowns. We have, discretizing the ODE the same way,
un+1 − 2un + un—
=f , un , ur ), n = 1, · · · , N − 1 (4.5)
1
(xn
h2 n

where ur = (u − u )/2h and ur = c and = d.

n n+1 n —1

0 ur N

The above formula not work for n = 0 or n = N because u—1 and uN+1 do not exist -
the stencil for the finite difference at point xn (using xn—1, xn, xn+1) would leave the domain.

Our last two equations come from the discretizing the boundary condition. It is not a
good idea to use a forward difference

y1 − y0
h ≈c
(why not?), and while a second-order forward difference could be used, it is not the ideal
solution. Instead, the typical approach is to use a centered difference to get

y1 − y—1
2h ≈ c
where y−1 denotes y(x—1) = y(−h). We now have an equation to extrapolate the value of y
18
at x—1, which yields an equation at the fictitious point x—1:

u—1 = u1 − 2hc. (4.6)

19
Similarly,
uN+1 = uN—1 + 2hc. (4.7)
Now we can apply the discretized formula (4.5) for the ODE at x0, since x−1 is now an
available point, leading to the single formula
un+1 − 2un + un—
1 =f , un , ur ), n = 0, · · · , N
(xn
h2 n

where u—1 and uN+1 are given by (4.6), (4.7). We now have a system of N + 1 equations for
the N + 1 unknowns that can be solved.

Essentially, the boundary conditions were used to pad the domain with ‘ghost points’ that
are not part of the solution, but are needed because the discretization extends past the end-
points. Depending on the problem, the trick here can be more delicate, but we will omit a
detailed discussion of the nuances here.

4.4 Midpoint discretization

Now consider the boundary value problem (in first-order form)

yr = F (x, y), x ∈ [0, b]

with m equations (so y : R → Rm) with linear boundary conditions

B0y(0) + Bby(b) = c.

Consider a grid {xn} with hn = xn+1 − xn and x0 = 0, xN = b and let xn+1/2 denote the
midpoint of [xn, xn+1]. We can use the midpoint formula to discretize the ODE - we use
a centered difference at each midpoint xn+1/2, which uses the values at xn and xn+1:
un+1 − un un+1 + un
=F , ), n = 0, · · · , N −
hn n+1/2 12
(x
noting that y at n + 1/2 has been approximated with an average.6

The boundary conditions provide one more set of equations, so there are mN equations
from the ODE and m from the BCs. This matches the m(N + 1) unknowns, so we have a
complete system that can be solved via the usual methods (e.g. Newton).

Note also that the resulting linear systems for Newton’s method will have a nice struc-
ture: except for the boundary conditions, the matrix will be banded - only elements a
certain distance from the main diagonal are non-zero.

If we let
∂F
Jn+1/2 =
(xn+1/2, un+1/2)
∂y
20
6
Note that there are other ways the discretization could be done; the version here is used as an example.

21
and
1 1
Vn = −I − hnJn+1/2, Wn = I − hnJn+1/2,
2 2
then the Jacobian for G(u) = 0 with

Gn(u) = un+1 − un − F (xn+1/2, un+1/2), n = 0, · · · , N − 1

and
GN (u) = B0u0 + BbuN − c
is given in block form by
 
V0 W0 0 ··· 0
0 V1 W1 ··· 0 
 .. .. ..
; 0 . . . .  .
0 ··· 0 VN—1 WN—1

B ‘ 
0 ··· 0 0 Bb 

If the boundary conditions are separated, i.e.

B˜0 y(0) B˜ b y(b) = 0, we can split the
=
boundary equations apart and put the x = 0 ones at the top rows and the x = b ones at the
bottom, preserving the band structure:
 ˜ 
B 0 0 ··· ··· 0
V0 W0 0 ··· 0
0 V1 W1 · · · 0
. . . .

; 0 .. .. .. .. 
0 ·· 0WN —1
VN —1

 · 

B ˜b
0 ·· 0 0
Because it is a ‘standard’ first order system, the method can be used on a variety of
problems (including the second-order BVP of the previous section). The primary advantage
is that the first order system allows us to discretize only first-order derivatives, which only
take two points (xn and xn+1). This small stencil makes using a non-uniform grid simple,
because the discretized ODE at each point xn+1/2 depends only only one of the grid
spacings (xn+1 − xn).

4.5 Upwind example

The ‘centered’ difference is not always the natural way to discretize. As an illustrative
example, consider
−λ 0
yr = y, y (0) = 1, y (3) = 1
1 2
0 λ
Observe that:

• The solution is the sum of two independent ‘modes’ ((y1(t), 0) + (0, y2(t)))
22
• Information propagates only to the right for y1 (from the boundary at x = 0)

23
 • Information propagates only to the left for y2 (from x = b)

The centered difference for the ODE at xn,

yn+1 − yn—1
ynr ≈
2h
depends on information in both directions (at indices n − 1 and n + 1). We might instead
want to have a method that has the same ‘propagation’ behavior as the solution. This means
that we should use backward differences for y1 and forward differences for y2, e.g.

(y )r
(y1)n − (y1)n—1 (y2)n+1 − (y2)n
≈
)r , (y ≈ .
1 n 2 n
h h
As seen in the stability discussion, this system is really just two IVPs in opposite directions.
If we use the upwind discretization, then we are using Backward Euler for both
components, so the method will be nice and stable even when λ is large.

If another discretization is used and λ is large, stiffness may be an issue.

The catch is that this approach is hard to generalize, since for a general system

yr = F (y)

it is difficult or impossible to separate into modes. This sort of discretization plays an

important role in solving PDEs with propagating waves (e.g. the wave equation, shocks
formed in gas flow, etc.), where there is a natural direction to the flow of information.

24
5 The linear algebra
We will not cover the numerical linear algebra here in detail 7. The focus here is on high-
lighting what one needs to know to use the algorithms.

A banded matrix A is a matrix such that

aij = 0 for − m1 ≤ j − i ≤ m2

where 0 < m1, m2 and m = m1 + m2 + 1 is the ‘bandwidth’, which should be much

smaller than the matrix dimensions. That is, A has non-zero diagonals only close to the
center diagonal. A common special case that we have seen is a tridiagonal matrix, for
which m1 = m − 2 = 1 and the bandwidth is three:

 γ1 0 ·· 0
 1β ·
α2 β2 γ2 .· .· 0 .
;  0 . .. . ..

. .
0 · · · αN —1 βN —1 γN —1

0 ··· 0 αN βN

An LU decomposition of a square matrix is a factorization

A = LU, L = lower triangular, U = upper triangular.

This factorization can be used to solve linear systems Ax = b in two steps:

1) (Factor) Compute L and U such that A = LU

2) (Solve) Solve Ly = b,then solve Ux = y (with forward/back substitution)

Once (1) is done, step (2) can be done for many RHS’s b without re-factoring the matrix,
so step (1) is often an initial ‘processing’ of the matrix to prepare it for solving linear systems.

In practice: When solving systems with LU decomposition, you should make use of a
specialized solver. This means, typically, calling code for (1) (e,g, L, U = lu(A)) and then
calling code for (2) to solve the system (e.g. x = solve(L,U,b)).

The LU factorization algorithm is straightforward (not described here), and the second step
is trivial. If A is N × N , then in general the factorization takes O(N 3) steps and the
second part takes O(N 2) steps. This means that in general, solving linear systems is expensive
and scales poorly as the size of the matrix increases.
See, for instance, Numerical Recipes for a good practical source or Golub and Van Loan’s Matrix com-
7

putations for details and theory

25
 5.1 LU for banded matrices
However, if A is banded then the solution process is much more efficient. The result:

(Banded LU) Let A be a banded matrix with m1 lower and m2 upper bands and bandwidth
m = m1 + m2 + 1. Then it can be factored as A = LU where

• L is lower triangular with m1 lower bands

• U is upper triangular with ≤ m1 + m2 upper bands

Computing A = LU takes O(m2N ) operations, as does solving Ax = b.

Because A only contains at most mN elements, we only need to store it in a ‘compressed’

form as an N × m array of diagonals. For example,
   
β 1 γ1 0 ··· 0 0 β 1 γ 1

α2 β2 γ2 ···  α.2 β.2 γ.2 

0 .
0 .. . .. . .. ..  =  . . . 
 γ  ⇒  
0 · · αN —1 βN —1 N —1 . . .
βN · αN β N 0
0 ·· 0 αN
Important note (padding): Because off-diagonals have less entries than the main
diago- nal, there will be some unused entries in the array. A choice must be made on how to
arrange the data (‘padding’ the array with zeros). The padded zeros are shown in red
above, using the convention that the j-th row of A is the j-th row of the compressed
array8.
As another example,

2 3 0 0 1 2 3
 4 5 6 7
1 4 5 6 7
 8 9 10
 
0 0 11 12

A=
  , m1 = 1, m2 = 2, =⇒ A (array) =  8 9 10 0
0   11

 12 0 0
Using the algorithms: A banded LU solver has two routines: a routine that computes
A = LU given A in compress form, and a routine that solves Ax = b given L, U (also
compressed) and b (see example code for typical usage). To solve Ax = b:

• Construct A in compressed form (be careful with indexing!)

• Call the routine to obtain L and U
• Call the routine to solve LUx = b

26
A brief note about other methods: If the matrix is sparse (has mostly non-zero en-
tries, e.g. an N × N matrix with O(N ) entries) but not banded, then other techniques must
8
Matlab uses the opposite convention for square matrices where columns line up; see the documentation
of spdiags for details.

27
be used. There are two common ones: iterative methods (e.g. GMRES) that construct
a sequence of approximations and do not require a specific pattern of non-zeros, and re-
ordering techniques, where rows and columns and other transforms are applied to move
the non-zeros near the diagonal.

Lastly, if the matrix is ‘almost banded’, such as the ‘periodic’ tridiagonal matrix

 γ1 0 ·· α1
β1 ·
α2 β2 γ
. .2 . .·
· 0
0 . . . 
 . .
 0 ·. · · α
 N —1 βN —1 γN —
 1 γN · · · 0 αN βN

That arises in solving problems with periodic boundary conditions, then the factorization
can be converted to a series of banded solves using the Sherman-Morrison formula, which
gives the inverse of a matrix A + uvT in terms of A—1 (where u, v are vectors and uvT is an
outer product).

28
6 More notes on BVPs
6.1 Quick example of a system
Let’s consider the BVP

(q(y)yr)r + cyr = f (y), y(0) = 0, y(1) = 3

with point xn = nh and h = 1/(N + 1). Suppose we want to convert to a first order system
and use the midpoint rule (this is useful e.g. if we have some non-uniform mesh to use).
For convenience, we use u and v as system variables, setting

u = y, v = qyr, u = (u, v), w = (u0, v0, u1, v1, · · · , uN+1, vN+1).

Remark (picking the right variable) It’s worth highlighting that the natural choice of
variables is (y, qyr) and not (y, yr). It’s often a good idea to choose your variables to simplify
the system; for physical problems, physical intuition can guide you. Here, qyr represents a
flux in a diffusion problem, and is a natural quantity to consider.

Then the problem to solve is

qur = v, vr = f (u) − cv, u(0) = u(1) = 0.

Applying the midpoint discretization (assuming a uniform step size for simplicity), we get
un+1 − un vn+1 − vn
qn+1/2 = , = — cvn+1/2 , n = 0, · · · , N
h h
vn+1/2 fn+1/2
along with u0 = 0, uN+1 = 3. This provides 2(N + 2) equations for all the values in w,
with two of them specified exactly. For the values at n + 1/2, an average is used. In some
cases, there are two choices, e.g.
1
qn+1/2 = (q(un) + q(un+1)) or qn+1/2 = q( un + un+1 )
2 2
both of which have the same order but slightly different properties (depending on the prob-
lem); not discussed here.

To be explicit, suppose q(y) = 1 and f (y) = ay and c = 0, so

yrr = ay, y(0) = 0, y(1) = 3

As before, we set u = y, v = yr and w = (u0, v0, · · · , uN+1, vN+1) and set

h
Gn(w) := un+1 − un − (vn + vn+1),
2 ah
Hn(w) := vn+1 − vn − (un + un+1), n = 0, · · · , N
2
29
along with u0 = uN+1 = 0. Stacking the equations in the order

u0 = 0, G0, H0, G1, H1, · · · GN , HN , · · · uN+1 = 0.

we get the following linear system:

   
1 0 0 ··· 0 u0 0
A0 B 0 0 ··· v0
 0   0 
0 A1 B 1 ·· ..   u   .
· 1 = .
.. . 
 . . .. ·· · .   .. 
. . .
 
   uN  
0 · · · 0 AN BN
 0
 3
+1

0 ··· 0 1 vN +1
0

30
where
 ∂Gn ∂Gn 
−1 −h/2
 ∂H ∂v 
An =  n ∂Hn  = — 2 −1
n ah

∂un ∂un ∂vn

 
∂Gn ∂Gn
1 −h
 ∂H
∂un+1 ∂v
∂H 
Bn =  n
n+1n
 = —ah
2
2
∂un+1 ∂vn+1
Note that to be efficient, one should ‘plug in’ the values of u0 and uN+1 and reduce the
system size by two; it is omitted here just to keep the structure simple.

6.2 Non-uniform meshes

(See Sec. 18.5 of Numerical Recipes, from which this section is adapted). We return to the
midpoint approach to the first order BVP

yr = F (x, y), ···.

With uniform spacing h , a midpoint discretization could be

un+1 − un un+1 + un
=F n+1/2 ,
h ).
(x
2
Often, it is better to have a non-uniform distribution of points where there are more points
in areas that need it (e.g. fast-varying parts of the solution needing more points).

One approach is to turn x into a dependent variable x(q), where the new independent
variable q is discretized with equal spacing. We choose the map from q → x.

Often we ∫ want a mesh density ρ(x) (a function that gives the density ∫of points on the
b x
grid, so ρ(x) = 1 and the number of points between x and x is ≈ N 2 ρ(x) dx).
a 1 2 x1

31
Example problem: Suppose, for the sake of an example, the density should be higher
where the slope of the first component y1 is large:
dy1
| | large =⇒ high density ρ(x). (6.1)
dx
Details: First, let us assume that q varies from 0 to 1 with uniform spacing ∆q. Since
dy dy dx dx
= =F
dq dx dq dq
we can discretize the above as an ODE with respect to q using the midpoint scheme. The
function x(q) could be chosen to get the density we want, but we may not know exactly
where the points should be in advance. Moreover, we need to ensure that q varies from 0 to
1 (our fixed interval). We seek q and a constant K so that
1 ∫ x
q(x) = ρ(ξ) dξ, q(1) = 1.
K 0

In addition, let ∫ x

Q(x) = ρ(ξ) dξ,

0
the anti-derivative of ρ. The value K (the normalization constant) is unknown, so we use
the constant to dependent variable trick. This gives ODEs (all with respect to q)
dQ dK dx K
= K, = =
dq 0, dq ρ
dq
which is then added to the main ODE
dy dx
=F
dq dq
leading to a system for (y, x, K, Q) with independent variable q. The density can then be
chosen with the desired features, for instance

ρ(x) = h0 dy1
+ C| |p
dx
where h0, C and p are constants to be configured. Note that ρ can depend on the solution
and its derivatives because it is part of the ODE, so it has access to y etc.

Obviously, a good choice of density will greatly improve efficiency (and often stability) by
using less points where they are not needed.

In addition, the strategy here is powerful because it provides a partly automated way to
select grid points, allowing one to use the solver to explore the features of the equation
without knowing everything in advance.
32
6.3 Galerkin
The Galerkin approach also represents the solution in terms of a basis, but imposes a very
different condition on the approximation. The idea of the Galerkan method, which is
typically what is known as ‘the’ finite element method, goes as follows. Again, for example,
let us consider the equation
yrr + cy = f (x), y(0) = y(1) = 0.
We consider approximation in terms of a basis {φk}M for an approximating subspace of the
k=0
set
{y : y(0) = y(1) = 0}.
The first step is to derive the weak form of the problem. We define the inner product
∫ 1
⟨f, g⟩ = 0 f (x)g(x)
dx
and take the inner product of the ODE with an arbitrary (smooth)9function ψ satisfying the
boundary conditions to get

— ⟨yr, ψr⟩ + c⟨y, ψ⟩ = ⟨f, ψ⟩ for all test functions ψ with ψ(0) = ψ(1) = 0 (6.2)

which is the weak form. Note that the BCs were used to integrate by parts and have otherwise
vanished. We now suppose our solution can be approximated in terms of the given basis,
M
Σ
y(x) ≈ akφk(x).
k=0

Obviously, we cannot require this approximation satisfies the weak form, since that must
hold for a continuous set of (arbitrary) ψ’s. Instead, we impose the condition that

(6.2) is satisfied for all test functions ψ in span{φk}

which leads to the equations

M M
Σ Σ
k n
−⟨ k=0akφr (x), φr ⟩ + c⟨k=0akφk(x), φn⟩ = ⟨f, φn⟩ for n = 0, · · · M.

Note that for (6.2) to hold for all ψ’s in the span of the φk’s, it suffices to require that it
holds for the basis functions.

In contrast to collocation, we do not require any condition to hold at a given mesh point; all
conditions are in terms of the weak form (involving only integrals over the interval).

This yields the system

V →a = f→
9
I am omitting some theory here that is not needed for the brief discussion of the numerics.

33
where
∫ 1
fk = ⟨f, φk⟩ = f (x)φk(x) dx,
0
Vnk = −⟨φr , φr ⟩ + c⟨φk, φn⟩.
k n

Now the subtle part: how do we choose the basis functions? A full discussion is beyond the
scope of the notes, but one simple choice should illustrate the point. Suppose we want an
approximation at a set of mesh points {xk} (for k = 0, · · · , M ).

We want basis functions with the property that

⟨φk, φn⟩ = 0 if k is not close to n

so that the resulting system will be sparse (in this case, banded). In addition, the basis
should be enough to represent any mesh function on the {xk}’s.

The piecewise linear basis we encountered in discussing splines is a good candidate. Re-
call these are piecewise linear ‘hats’
(
(line from (xk—1, 0) to ((xk, 1)) xk—1 < x < xk
φk (x) = (line from (xk, 1) to ((xk+1, 0)) xk < x < xk+1

for k = 1, · · · , M − 1. This is the basis we would use, as it spans mesh functions on the given
points that satisfy the boundary conditions. We can then compute the quantities ⟨φr , φr ⟩
k n
and so on and solve the system.

The method can be generalized by taking other bases of ‘elements’ - overlapping functions
with small support - to obtain better accuracy, stability or other properties. The method
is popular for complicated problems in multiple dimensions because it does not care much
about the location of the mesh points, allowing one to write general code that can handle
complicated geometries.

34