Numerical Methods for Differential Equations

Chapter 1: Initial value problems in ODEs

Gustaf Soderlind and Carmen Arevalo

Numerical Analysis, Lund University

Textbooks: A First Course in the Numerical Analysis of Differential Equations, by Arieh Iserles
and Introduction to Mathematical Modelling with Differential Equations, by Lennart Edsberg

c Gustaf Soderlind, Numerical Analysis, Mathematical Sciences, Lund University, 2008-09

Numerical Methods for Differential Equations p. 1/52

 Chapter 1: contents

Course contents
Introduction to initial value problems
The explicit Euler method
Convergence
Order of consistency
The trapezoidal rule
Theta methods
Numerical tests
The linear test equation and numerical stability
Stiff equations
Numerical Methods for Differential Equations p. 2/52
 What is Numerical Analysis?

Categories of mathematical problems

Category Algebra Analysis

linear computable not computable
nonlinear not computable not computable

Algebra: Only finite constructions

Analysis: Limits, derivatives, integrals etc. (transfinite)
Computable: the exact solution can be obtained in a
finite number of operations

Numerical Methods for Differential Equations p. 3/52

 What is Numerical Analysis?

Category Algebra Analysis

linear computable not computable
nonlinear not computable not computable

With numerical methods, problems from all four categories

can be solved:

Numerical analysis aims to construct and analyze

quantitative methods for the automatic computation of
approximate solutions to mathematical problems.

Goal: Construction of mathematical software

Numerical Methods for Differential Equations p. 4/52
 0. What will we study in this course?

To solve a differential equation analytically we look for a

differentiable function that satisfies the equation
Large, complex and nonlinear systems cannot be solved
analytically

Numerical Methods for Differential Equations p. 5/52

 0. What will we study in this course?

To solve a differential equation analytically we look for a

differentiable function that satisfies the equation
Large, complex and nonlinear systems cannot be solved
analytically

Instead, we compute numerical solutions with standard

methods and software

To solve a differential equation numerically we generate a

sequence {yk }N k=0 of pointwise approximations to the
analytical solution: y(tk ) yk

Numerical Methods for Differential Equations p. 5/52

 Some differential equations we will solve

Initial value problems (IVP) first-order equations;

higher-order equations; systems of differential
equations

Numerical Methods for Differential Equations p. 6/52

 Some differential equations we will solve

Initial value problems (IVP) first-order equations;

higher-order equations; systems of differential
equations

Boundary value problems (BVP) two-point boundary

value problems; Sturm-Liouville eigenvalue problems

Numerical Methods for Differential Equations p. 6/52

 Some differential equations we will solve

Initial value problems (IVP) first-order equations;

higher-order equations; systems of differential
equations

Boundary value problems (BVP) two-point boundary

value problems; Sturm-Liouville eigenvalue problems

Partial differential equations (PDE) the diffusion

equation; the advection equation; the wave equation

Numerical Methods for Differential Equations p. 6/52

 Some differential equations we will solve

Initial value problems (IVP) first-order equations;

higher-order equations; systems of differential
equations

Boundary value problems (BVP) two-point boundary

value problems; Sturm-Liouville eigenvalue problems

Partial differential equations (PDE) the diffusion

equation; the advection equation; the wave equation

Applications in all three areas

Numerical Methods for Differential Equations p. 6/52

 Initial value problems: examples

A first-order equation: a simple equation without a known

analytical solution
dy t2
=ye , y(0) = y0
dt

Numerical Methods for Differential Equations p. 7/52

 Initial value problems: examples

A first-order equation: a simple equation without a known

analytical solution

A second-order equation: motion of a pendulum

g
(t) + sin (t) = 0, (0) = 0 , (0) = 0
L
where (t) is the angle, g is the gravitational constant and
L is the pendulum length

Numerical Methods for Differential Equations p. 7/52

 Initial value problems: examples

A first-order equation: a simple equation without a known

analytical solution

A second-order equation: pendulum equation

A system of equations: the predatorprey model

y1 (t) = k1 y1 (t) k2 y1 (t) y2 (t)
y2 (t) = k3 y1 (t) y2 (t) k4 y2 (t)

where y1 (t) is the population of the prey species and y2 (t)

is the population of the predator species at time t

Numerical Methods for Differential Equations p. 7/52

 Boundary value problems: examples

Second-order two-point BVP: the electrostatic potential

u(r) between two concentric metal spheres satisfies

d2 u 2 du
2
+ = 0, u(R1 ) = V1 , u(R2 ) = 0
dr r dr
at the distance r from the center; R1 and R2 are the radii of
the two spheres

Numerical Methods for Differential Equations p. 8/52

 Boundary value problems: examples

Second-order two-point BVP: the electrostatic potential

u(r) between two concentric metal spheres satisfies

d2 u 2 du
2
+ = 0, u(R1 ) = V1 , u(R2 ) = 0
dr r dr
at the distance r from the center; R1 and R2 are the radii of
the two spheres

A Sturm-Liouville eigenproblem: Euler buckling of column

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

Find eigenvalues and eigenfunctions y

Numerical Methods for Differential Equations p. 8/52

 Partial differential equations: examples

The heat equation

ut (x, t) = uxx (x, t), x [0, a), t (0, b)
u(x, 0) = f (x), x [0, a]
u(0, t) = c1 , u(a, t) = c2 , t [0, b]

is a parabolic PDE modelling e.g. the temperature in an

insulated rod with constant temperatures c1 and c2 at its
ends, and initial temperature distribution f (x)

Numerical Methods for Differential Equations p. 9/52

 Partial differential equations: examples

The heat equation

The wave equation

utt (x, t) = uxx (x, t), x (0, a), t (0, b)
u(0, t) = 0, u(a, t) = 0 for t [0, b]
u(x, 0) = f (x) for x [0, a]
ut (x, 0) = g(x) for x (0, a)

is a hyperbolic PDE e.g. modelling the displacement u of a

vibrating elastic string fixed at x = 0 and x = a

Numerical Methods for Differential Equations p. 9/52

 Some applications in ODEs

Initial value problems

mechanics M q = F (q)
electrical circuits C v = I(v)
chemical reactions c = f (c)

Boundary value problems

materials u = M (u)/EI
microphysics ~
2m = E
eigenmodes u + u = 0

Numerical Methods for Differential Equations p. 10/52

 1. Initial value problems

Standard formulation

y = f (t, y) ; y(0) = y0

Scalar equation: y R; f scalar;

Systems of ODEs: y Rm ; f vector-valued.

Theorem (Cauchy & Peano) A solution to y = f (t, y)

exists if f is continuous.

Numerical Methods for Differential Equations p. 11/52

 Existence
Note: Continuity is not sufficient to guarantee uniqueness!

Example:

y = 2 y; y(0) = 0 with solution
y(t) = 0, t and y(t) = (t )2 , t >

This may happen in real physical models!

Consider the energy of a falling particle of mass m:

1
my 2 + mgy = E := 0
2

This yields y = 2g y with y(0) = 0
Numerical Methods for Differential Equations p. 12/52
 Existence and uniqueness

Theorem If f (t, y) is continuous on [0, T ] and satisfies the

Lipschitz condition

kf (t, u) f (t, v)k L[f ] ku vk

for all u, v and some Lipschitz constant L[f ] < , there

exists a unique solution to the initial value problem

y = f (t, y)

on [0, T ] for every initial value y(0) = y0

Numerical Methods for Differential Equations p. 13/52

 Existence and uniqueness . . .

Note Most problems do not satisfy a Lipschitz condition

on all of Rm

Example The problem

y = y 2 ; y(0) = y0 > 0
has solution
y0
y(t) = 1y0 t

The solution blows up at t = 1/y0 (Finite escape time)

Numerical Methods for Differential Equations p. 14/52

 Existence and uniqueness . . .

Other problems always satisfy Lipschitz conditions on Rm

Example
y = Ay; y(0) = y0

Lipschitz constant:
kAuAvk kAyk
L[f ] = maxu6=v kuvk
= maxy6=0 kyk
= kAk

The matrix norm kAk is a Lipschitz constant for f (y) = Ay

Numerical Methods for Differential Equations p. 15/52

 Standard form: x = F (t, x)
Example

y = f (t, y, y ); y(0) = y0 ; y (0) = y0

Standard substitution Introduce new variables

x1 = y
x2 = y
Then we get a system of first order equations
x1 = x2
x2 = f (t, x1 , x2 )
with x1 (0) = y0 and x2 (0) = y0

Numerical Methods for Differential Equations p. 16/52

 2. The Explicit Euler method (1768)
y = f (t, y); y(t0 ) = y0

Replace derivative by finite difference approximation

y(tn + h) y(tn )
y (tn )
h
Let {un } denote the numerical approximation to {y(tn )} and
h = tn+1 tn denote the time-step.

Numerical Methods for Differential Equations p. 17/52

 2. The Explicit Euler method (1768)
y = f (t, y); y(t0 ) = y0

Replace derivative by finite difference approximation

y(tn + h) y(tn )
y (tn )
h
Let {un } denote the numerical approximation to {y(tn )} and
h = tn+1 tn denote the time-step. Compute un from
un+1 un
= f (tn , un ), u0 = y0
h

Numerical Methods for Differential Equations p. 17/52

 2. The Explicit Euler method (1768)
y = f (t, y); y(t0 ) = y0

Replace derivative by finite difference approximation

y(tn + h) y(tn )
y (tn )
h
Let {un } denote the numerical approximation to {y(tn )} and
h = tn+1 tn denote the time-step. Compute un from
un+1 un
= f (tn , un ), u0 = y0
h
The Explicit Euler method Given u0 , t0 and h, compute
un+1 = un + hf (tn , un )
tn+1 = tn + h
Numerical Methods for Differential Equations p. 17/52
 Explicit Euler: Taylor series expansion

y (t) = f (t, y), y(t0 ) = y0

Use Taylor series expansion
2
h
y(t + h) = y(t) + hy (t) + y ()
2!
= y(t) + hf (t, y(t)) + O(h2 )
y(t + h) y(t) + hf (t, y(t))

Numerical Methods for Differential Equations p. 18/52

 Explicit Euler: Taylor series expansion

y (t) = f (t, y), y(t0 ) = y0

Use Taylor series expansion
2
h
y(t + h) = y(t) + hy (t) + y ()
2!
= y(t) + hf (t, y(t)) + O(h2 )
y(t + h) y(t) + hf (t, y(t))

Construct the numerical method (drop higher order terms)

un+1 = un + hf (tn , un ); u0 = y0
tn+1 = tn + h

Numerical Methods for Differential Equations p. 18/52

 Explicit Euler: Taylor series expansion

y (t) = f (t, y), y(t0 ) = y0

Use Taylor series expansion
2
h
y(t + h) = y(t) + hy (t) + y ()
2!
= y(t) + hf (t, y(t)) + O(h2 )
y(t + h) y(t) + hf (t, y(t))

Construct the numerical method (drop higher order terms)

un+1 = un + hf (tn , un ); u0 = y0
tn+1 = tn + h

Explicit Euler!
Numerical Methods for Differential Equations p. 18/52
 Explicit Euler: graphic interpretation
Take a step of size h in the direction of the tangent
1.2

1.15

1.1

1.05

0.95
y

0.9

0.85

0.8

0.75

0.7
0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
t

Each step introduces an error and ends up on a different

solution trajectory (dashed curves)
Numerical Methods for Differential Equations p. 19/52
 Explicit Euler: an example
Compute numerical solution to y = y cos t; y(0) = 1; t [0, 2]

Solution Choose a stepsize h = 2/N ; N = 24 h = /12

Recursion un+1 = un + h (un cos tn ); u0 = 1

tn+1 = tn + h; t0 = 0
1.8 1.8

1.6 1.6

1.4 1.4

1.2 1.2

1 1

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2
0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7

The numerical solution is a sequence of points (tn , un )

Numerical Methods for Differential Equations p. 20/52
 3. Convergence
Analytical and numerical solutions for h = /8 and /128
3 3

N=64 N=64*16

2.5 2.5

2 2

1.5 1.5

1 1

0.5 0.5

0 0
0 5 10 15 20 25 0 5 10 15 20 25

As h 0 the numerical solution approaches the exact solution

Numerical Methods for Differential Equations p. 21/52

 3. Convergence
Analytical and numerical solutions for h = /8 and /128
3 3

N=64 N=64*16

2.5 2.5

2 2

1.5 1.5

1 1

0.5 0.5

0 0
0 5 10 15 20 25 0 5 10 15 20 25

As h 0 the numerical solution approaches the exact solution

A method is convergent if, for every ODE with a Lipschitz
function f and every fixed T , with T = N h, it holds that

lim kyN,h y(T )k = 0

N

Numerical Methods for Differential Equations p. 21/52

 Local and global errors
1.15

yn+1

1.1

yn

1.05
y

yn+1

1

y(tn ) y(tn+1 ) y(t)
0.95
1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7 1.75 1.8
t

Global error en = yn y(tn ) and en+1 = yn+1 y(tn+1 )

Local error ln+1 = yn+1 y(tn+1 )

Numerical Methods for Differential Equations p. 22/52

 The local error of the explicit Euler method
Insert exact data

yn+1 = y(tn ) + hf (tn , y(tn ))

Local error definition ln+1 = yn+1 y(tn+1 ) implies residual

y(tn+1 ) = y(tn ) + hf (tn , y(tn )) ln+1

h2
Taylor series y(tn+1 ) = y(tn ) + hy (tn ) + 2
y (tn ) + ...

2
Local error for explicit Euler ln+1 h2 y (tn )
(along exact solution, solid curve)
Numerical Methods for Differential Equations p. 23/52
 Error propagation

Explicit Euler (numerical solution)

yn+1 = yn + hf (tn , yn )

Subtract Taylor series expansion of exact solution

h2
y(tn+1 ) = y(tn ) + hf (tn , y(tn )) + y (tn ) + . . .
2
Global error recursion

en+1 = en + hf (tn , y(tn ) + en ) hf (tn , y(tn )) + ln+1

Numerical Methods for Differential Equations p. 24/52

 Error propagation . . .

en+1 = en + hf (tn , y(tn ) + en ) hf (tn , y(tn )) + ln+1

Take norms and use Lipschitz condition:

ken+1 k ken k + hL[f ] ken k + kln+1 k

Lemma If {an }, a0 = 0, is a sequence of non-negative

numbers satisfying

an+1 (1 + h )an + ch2 for > 0, then

c
an h [(1 + h )n 1], n = 0, 1, . . .

Numerical Methods for Differential Equations p. 25/52

 Convergence of Eulers method: Theorem

Theorem The explicit Euler method is convergent

Numerical Methods for Differential Equations p. 26/52

 Convergence of Eulers method: Theorem

Theorem The explicit Euler method is convergent

Proof Suppose f sufficiently differentiable. Given h and a

fixed T = N h, let en,h = yn,h y(tn )

Numerical Methods for Differential Equations p. 26/52

 Convergence of Eulers method: Theorem

Theorem The explicit Euler method is convergent

Proof Suppose f sufficiently differentiable. Given h and a

fixed T = N h, let en,h = yn,h y(tn )

Apply the lemma to global error recursion, to get

c
ken,h k h [(1 + h L[f ])n 1], n = 0, 1, . . .
L[f ]

with
c = max kln k/h2 max ky k/2
n t

Numerical Methods for Differential Equations p. 26/52

 Convergence of Eulers method . . .

As (1 + h L[f ])n < enhL[f ] eT L[f ] , we have for n N ,

c
ken,h k h (eT L[f ] 1)
L[f ]

So
ken,h k C(T ) h
implies convergence because

lim ken,h k = 0
h0

Note 1) The global error can be made arbitrarily small!

2) The error is way too large for practical purposes!
3) Better bounds can be obtained
Numerical Methods for Differential Equations p. 27/52
 Theoretical error bound
Example
y = 100y, y(0) = 1.
Then L[f ] = 100 and the exact solution is y(t) = e100t with
y (t) = 1002 e100t , so c = 1002 /2, with bound

1002
ken,h k h (e100T 1)
2 100

Numerical Methods for Differential Equations p. 28/52

 Theoretical error bound
Example
y = 100y, y(0) = 1.
Then L[f ] = 100 and the exact solution is y(t) = e100t with
y (t) = 1002 e100t , so c = 1002 /2, with bound

1002
ken,h k h (e100T 1)
2 100

Error estimate at T = 1 is ken,h k 50 h e100 1.4 1045 h !

Numerical Methods for Differential Equations p. 28/52

 Theoretical error bound
Example
y = 100y, y(0) = 1.
Then L[f ] = 100 and the exact solution is y(t) = e100t with
y (t) = 1002 e100t , so c = 1002 /2, with bound

1002
ken,h k h (e100T 1)
2 100

Error estimate at T = 1 is ken,h k 50 h e100 1.4 1045 h !

But yn = (1 100h)n , so for h < 1/50, at T = 1 ,

Actual error is ken,h k = |(1 100/N )N e100 | 3.7 1044 h !

Numerical Methods for Differential Equations p. 28/52

 Computational test y = (y sin t) + cos t

= 0.2 , with initial condition y(/4) = 1/ 2
1.2 1.2

1.1 1.1

1 1

0.9 0.9
y

y
0.8 0.8

0.7 0.7

0.6 0.6

0.5 0.5
0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
t t

h = /10 h = /20

Note Local error O(h2 ), global error O(h)!

Numerical Methods for Differential Equations p. 29/52
 4. Order of consistency
Given a method in the form yn+1 = n (f, h, y0 , y1 , . . . , yn )
Insert exact data y(t0 ), y(t1 ), . . . , y(tn )

Definition The order of consistency is p if

y(tn+1 ) n (f, h, y(t0 ), y(t1 ), . . . , y(tn )) = O(hp+1 ), n

as h 0, for every analytic f

Theorem The local error is then O(hp+1 )

Numerical Methods for Differential Equations p. 30/52

 4. Order of consistency
Given a method in the form yn+1 = n (f, h, y0 , y1 , . . . , yn )
Insert exact data y(t0 ), y(t1 ), . . . , y(tn )

Definition The order of consistency is p if

y(tn+1 ) n (f, h, y(t0 ), y(t1 ), . . . , y(tn )) = O(hp+1 ), n

as h 0, for every analytic f

Theorem The local error is then O(hp+1 )

Alternative The order of consistency is p if the formula is

exact for all polynomials y = P (t) of degree p or less

Numerical Methods for Differential Equations p. 30/52

 Order of consistency of Eulers method
Example Explicit Euler n (f, h, y0 , . . . , yn ) = yn + h f (tn , yn )

Numerical Methods for Differential Equations p. 31/52

 Order of consistency of Eulers method
Example Explicit Euler n (f, h, y0 , . . . , yn ) = yn + h f (tn , yn )
Expanding in Taylor series,
y(tn+1 ) [y(tn ) + h f (tn , y(tn ))] = O(h2 )

so the methods consistency order is one

Numerical Methods for Differential Equations p. 31/52

 Order of consistency of Eulers method
Example Explicit Euler n (f, h, y0 , . . . , yn ) = yn + h f (tn , yn )
Expanding in Taylor series,
y(tn+1 ) [y(tn ) + h f (tn , y(tn ))] = O(h2 )

so the methods consistency order is one

Alternatively, suppose y(t) = 1 with f = y = 0. Then

1 = y(tn+1 ) = y(tn ) + hf (tn , y(tn )) = 1 + h 0 = 1 exact!

Numerical Methods for Differential Equations p. 31/52

 Order of consistency of Eulers method
Example Explicit Euler n (f, h, y0 , . . . , yn ) = yn + h f (tn , yn )
Expanding in Taylor series,
y(tn+1 ) [y(tn ) + h f (tn , y(tn ))] = O(h2 )

so the methods consistency order is one

Alternatively, suppose y(t) = 1 with f = y = 0. Then

1 = y(tn+1 ) = y(tn ) + hf (tn , y(tn )) = 1 + h 0 = 1 exact!

Next take 1st degree polynomial y(t) = t with f = y = 1, then

h = y(tn+1 ) = y(tn ) + hf (tn , y(tn )) = 0 + h 1 = h exact!

Numerical Methods for Differential Equations p. 31/52

 Order of consistency of Eulers method
Example Explicit Euler n (f, h, y0 , . . . , yn ) = yn + h f (tn , yn )
Expanding in Taylor series,
y(tn+1 ) [y(tn ) + h f (tn , y(tn ))] = O(h2 )

so the methods consistency order is one

Alternatively, suppose y(t) = 1 with f = y = 0. Then

1 = y(tn+1 ) = y(tn ) + hf (tn , y(tn )) = 1 + h 0 = 1 exact!

Next take 1st degree polynomial y(t) = t with f = y = 1, then

h = y(tn+1 ) = y(tn ) + hf (tn , y(tn )) = 0 + h 1 = h exact!

For a 2nd degree polynomial, y(t) = t2 with f = y = 2t. Then

h2 = y(tn+1 ) = y(tn ) + hf (tn , y(tn )) = 0 + h 0 = 0 6= h2 p = 1 !

Numerical Methods for Differential Equations p. 31/52

 5. The trapezoidal rule
Explicit Euler linearizes y at tn with a slope of y (tn )

Numerical Methods for Differential Equations p. 32/52

 5. The trapezoidal rule
Explicit Euler linearizes y at tn with a slope of y (tn )
Instead, linearize with a slope equal to the average of y (tn )
and y (tn+1 )
1
y(t) y(tn ) + (t tn ) [f (tn , y(tn )) + f (tn+1 , y(tn+1 ))]
2

Numerical Methods for Differential Equations p. 32/52

 5. The trapezoidal rule
Explicit Euler linearizes y at tn with a slope of y (tn )
Instead, linearize with a slope equal to the average of y (tn )
and y (tn+1 )
1
y(t) y(tn ) + (t tn ) [f (tn , y(tn )) + f (tn+1 , y(tn+1 ))]
2

The numerical method is the trapezoidal rule

tn+1 = tn + h
h
yn+1 = yn + [f (tn , yn ) + f (tn+1 , yn+1 )]
2

Numerical Methods for Differential Equations p. 32/52

 5. The trapezoidal rule
Explicit Euler linearizes y at tn with a slope of y (tn )
Instead, linearize with a slope equal to the average of y (tn )
and y (tn+1 )
1
y(t) y(tn ) + (t tn ) [f (tn , y(tn )) + f (tn+1 , y(tn+1 ))]
2

The numerical method is the trapezoidal rule

tn+1 = tn + h
h
yn+1 = yn + [f (tn , yn ) + f (tn+1 , yn+1 )]
2
The method is implicit. Nonlinear equation solving required

Numerical Methods for Differential Equations p. 32/52

 Trapezoidal rule: Order and convergence
Order of consistency (sketch):

h
y(tn+1 ) {y(tn ) + [f (tn , y(tn )) + f (tn+1 , y(tn+1 ))]}
2
2
h
= y(tn ) + h y (tn ) + y (tn ) + O(h3 )
2
h
{y(tn ) + (y (tn ) + [y (tn ) + h y (tn ) + O(h2 )])}
2
= O(h3 )

Numerical Methods for Differential Equations p. 33/52

 Trapezoidal rule: Order and convergence
Order of consistency (sketch):

h
y(tn+1 ) {y(tn ) + [f (tn , y(tn )) + f (tn+1 , y(tn+1 ))]}
2
2
h
= y(tn ) + h y (tn ) + y (tn ) + O(h3 )
2
h
{y(tn ) + (y (tn ) + [y (tn ) + h y (tn ) + O(h2 )])}
2
= O(h3 )

The method is of order two

Numerical Methods for Differential Equations p. 33/52

 Trapezoidal rule: Order and convergence
Order of consistency (sketch):

h
y(tn+1 ) {y(tn ) + [f (tn , y(tn )) + f (tn+1 , y(tn+1 ))]}
2
2
h
= y(tn ) + h y (tn ) + y (tn ) + O(h3 )
2
h
{y(tn ) + (y (tn ) + [y (tn ) + h y (tn ) + O(h2 )])}
2
= O(h3 )

The method is of order two

Theorem The trapezoidal rule is convergent

(No proof given here)
Numerical Methods for Differential Equations p. 33/52
Trapezoidal rule: the dramatic effect of 2nd order
Compute numerical solution to y = y cos t; y(0) = 1; t [0, 8].
Choose h = /12. Note that f is linear in y. The recursion is
h
1 2 cos tn
un+1 = h
un ; u0 = 1
1+ 2 cos tn+1
3 3
N=96 Trapezoidal rule N=96 Explicit Euler

2.5 2.5

2 2

1.5 1.5

1 1

0.5 0.5

0 0
0 5 10 15 20 25 0 5 10 15 20 25

Solutions with Trapezoidal rule and Explicit Euler

Numerical Methods for Differential Equations p. 34/52
 The implicit midpoint rule

We can also approximate the derivative y (t) by taking the

average of tn and tn+1 as well as yn and yn+1 :

tn + tn+1 yn + yn+1
y (t) f ( , ), t [tn , tn+1 ]
2 2

Numerical Methods for Differential Equations p. 35/52

 The implicit midpoint rule

We can also approximate the derivative y (t) by taking the

average of tn and tn+1 as well as yn and yn+1 :

tn + tn+1 yn + yn+1
y (t) f ( , ), t [tn , tn+1 ]
2 2

The resulting method is the 2nd order implicit midpoint method

tn+1 = tn + h
tn + tn+1 yn + yn+1
yn+1 = yn + h f ( , )
2 2

Numerical Methods for Differential Equations p. 35/52

 6. Theta methods

We construct methods that linearize y at tn with a slope equal to

a convex combination of y (tn ) and y (tn+1 ):

yn+1 = yn + h [f (tn , yn ) + (1 )f (tn+1 , yn+1 )], [0, 1]

Explicit Euler: = 1
Trapezoidal rule (implicit): = 1/2
Implicit Euler: = 0 yn+1 = yn + h f (tn+1 , yn+1 )

Numerical Methods for Differential Equations p. 36/52

 Theta methods: Order and convergence

Use Taylor expansion to get

y(tn+1 ) y(tn ) h [f (tn , y(tn )) + (1 )f (tn+1 , y(tn+1 ))]
1 1 2
= ( )h2 y (tn ) + ( )h3 y (tn ) + O(h4 )
2 2 3

Numerical Methods for Differential Equations p. 37/52

 Theta methods: Order and convergence

Use Taylor expansion to get

y(tn+1 ) y(tn ) h [f (tn , y(tn )) + (1 )f (tn+1 , y(tn+1 ))]
1 1 2
= ( )h2 y (tn ) + ( )h3 y (tn ) + O(h4 )
2 2 3
If = 1/2 the method is of order 2; otherwise it is of order 1

Numerical Methods for Differential Equations p. 37/52

 Theta methods: Order and convergence

Use Taylor expansion to get

y(tn+1 ) y(tn ) h [f (tn , y(tn )) + (1 )f (tn+1 , y(tn+1 ))]
1 1 2
= ( )h2 y (tn ) + ( )h3 y (tn ) + O(h4 )
2 2 3
If = 1/2 the method is of order 2; otherwise it is of order 1

Theorem (without proof) The methods are convergent

Numerical Methods for Differential Equations p. 37/52

 Implicit methods

Implicit Euler yn+1 = yn + h f (tn+1 , yn+1 )

We need to solve a nonlinear equation to compute yn+1
The extra cost is motivated if we can take larger steps

There are some problems where implicit methods can take

enormous time steps without losing accuracy!

We will return to how to solve nonlinear equations

Numerical Methods for Differential Equations p. 38/52

 7. Tests: Explicit Euler y = (y sin t) + cos t

= 0.2 , with initial condition y(/4) = 1/ 2
1.2 1.2

1.1 1.1

1 1

0.9 0.9
y

y
0.8 0.8

0.7 0.7

0.6 0.6

0.5 0.5
0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
t t

h = /10 h = /20

Local error O(h2 ), global error O(h)!

Numerical Methods for Differential Equations p. 39/52
 Tests: Implicit Euler y = (y sin t) + cos t

= 0.2 , with initial condition y(/4) = 1/ 2
1.2 1.2

1.1 1.1

1 1

0.9 0.9
y

y
0.8 0.8

0.7 0.7

0.6 0.6

0.5 0.5
0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
t t

h = /10 h = /20

Local error O(h2 ), global error O(h)!

Numerical Methods for Differential Equations p. 40/52
 Tests: Implicit Euler y = (y sin t) + cos t

= 10 , with initial condition y(/4) = 1/ 2
1.2 1.2

1.1 1.1

1 1

0.9 0.9
y

y
0.8 0.8

0.7 0.7

0.6 0.6

0.5 0.5
0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
t t

h = /10 h = /20

Local error = O(h2 ), global error = O(h), but difficult to see!

Numerical Methods for Differential Equations p. 41/52
 Tests: Explicit Euler y = (y sin t) + cos t

= 10 , with initial condition y(/4) = 1/ 2
1.2

1.1

0.9

y
0.8

0.7

0.6

0.5
0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
t

h = /20

Numerical Methods for Differential Equations p. 42/52

 Tests: Explicit Euler y = (y sin t) + cos t

= 10 , with initial condition y(/4) = 1/ 2
1.2 1.2

1.1 1.1

1 1

0.9 0.9
y

y
0.8 0.8

0.7 0.7

0.6 0.6

0.5 0.5
0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
t t

h = /10 h = /20

NUMERICAL INSTABILITY! Stability for h small

Numerical Methods for Differential Equations p. 42/52
 8. The linear test equation

Defintion The linear test equation is

y = y; y(0) = 1, t 0, where C

Numerical Methods for Differential Equations p. 43/52

 8. The linear test equation

Defintion The linear test equation is

y = y; y(0) = 1, t 0, where C

As y(t) = et , we have

|y(t)| K Re() 0

Mathematical stability Bounded solutions if Re() 0

Numerical Methods for Differential Equations p. 43/52

 8. The linear test equation

Defintion The linear test equation is

y = y; y(0) = 1, t 0, where C

As y(t) = et , we have

|y(t)| K Re() 0

Mathematical stability Bounded solutions if Re() 0

When does a numerical method have the same property?

Does Re() 0 imply numerical stability?

Numerical Methods for Differential Equations p. 43/52
 Numerical stability: The stability region

Definition The stability region D of a method is the set of all

h C such that |yn | K when the method is applied to the test
equation

Numerical Methods for Differential Equations p. 44/52

 Numerical stability: The stability region

Definition The stability region D of a method is the set of all

h C such that |yn | K when the method is applied to the test
equation

Example For Eulers method, yn+1 = (1 + h)yn , so yn remains

bounded if and only if |1 + h| 1

Numerical Methods for Differential Equations p. 44/52

 Numerical stability: The stability region

Definition The stability region D of a method is the set of all

h C such that |yn | K when the method is applied to the test
equation

Example For Eulers method, yn+1 = (1 + h)yn , so yn remains

bounded if and only if |1 + h| 1
4

3
DEuler = {z C : |1 + z| 1}
2

4
4 3 2 1 0 1 2 3 4

Numerical Methods for Differential Equations p. 44/52

 A-stability
Definition A method is called A-stable if its stability region
contains C , i.e. C {z C : Re(z) 0} D.

Numerical Methods for Differential Equations p. 45/52

 A-stability
Definition A method is called A-stable if its stability region
contains C , i.e. C {z C : Re(z) 0} D.

For the trapezoidal rule

(   )
1 + 1z 
2 
DTR = zC: 1 C

 1 12 z 


so if Re(z) 0 then yn remains bounded for any h > 0

Numerical Methods for Differential Equations p. 45/52

 A-stability
Definition A method is called A-stable if its stability region
contains C , i.e. C {z C : Re(z) 0} D.

For the trapezoidal rule

(   )
1 + 1z 
2 
DTR = zC: 1 C

 1 12 z 


so if Re(z) 0 then yn remains bounded for any h > 0

The explicit Euler method is not A-stable but the implicit Euler
method and the trapezoidal rule are A-stable

Usually: If the original problem is stable, then an A-stable

method will replicate that behavior numerically
Numerical Methods for Differential Equations p. 45/52
 Relevance of the linear test equation

When we have a system of equations

y = Ay

and A has a full set of eigenvectors with eigenvalues

= diag(1 , 2 , . . . , N ), then
A can be diagonalized, T 1 AT =
Putting y = T x transforms the system to x = x (scalar
systems)
. . . equal to xk = k xk for k = 1, . . . , d
The exact solution y(t) = etA y0 is stable if and only if
Re(k ) 0 for all k = 1, . . . , d

Numerical Methods for Differential Equations p. 46/52

 Relevance of the linear test equation . . .
Applying a method to y = Ay gives (say)

yn+1 = (I + hA)yn

Numerical Methods for Differential Equations p. 47/52

 Relevance of the linear test equation . . .
Applying a method to y = Ay gives (say)

yn+1 = (I + hA)yn

Noting that T 1 (I + hA)T = I + h, putting yn = T xn produces

xn+1 = (I + h)xn

Numerical Methods for Differential Equations p. 47/52

 Relevance of the linear test equation . . .
Applying a method to y = Ay gives (say)

yn+1 = (I + hA)yn

Noting that T 1 (I + hA)T = I + h, putting yn = T xn produces

xn+1 = (I + h)xn

The same as applying the method to the diagonalized system

x = x

Numerical Methods for Differential Equations p. 47/52

 Relevance of the linear test equation . . .
Applying a method to y = Ay gives (say)

yn+1 = (I + hA)yn

Noting that T 1 (I + hA)T = I + h, putting yn = T xn produces

xn+1 = (I + h)xn

The same as applying the method to the diagonalized system

x = x

Diagonalization and discretization commute! Stability condition

from linear test equation applies to all diagonalizable systems!
Numerical Methods for Differential Equations p. 47/52
 9. Stiff ODEs

Example Solve y = (y sin t) + cos t with = 50

Solution
Particular: yP (t) = sin t
Homogeneous: yH (t) = et
General: y(t) = e(tt0 ) (y(t0 ) sin t0 ) + sin t

Study the flow of this equation and numerical solutions

Numerical Methods for Differential Equations p. 48/52

 Flow (solution trajectories)

1.2

1.1

0.9

0.8
y

0.7

0.6

0.5

0.4
0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4
t

Numerical Methods for Differential Equations p. 49/52

 Stiffness. . .

With Explicit Euler the solution approaches sin t as t ,

1
as the exact solution does, only if 0 < h <
25

Numerical Methods for Differential Equations p. 50/52

 Stiffness. . .

With Explicit Euler the solution approaches sin t as t ,

1
as the exact solution does, only if 0 < h <
25
This means h must be kept small, not to keep errors small,
but to have a stable numerical solution

Numerical Methods for Differential Equations p. 50/52

 Stiffness. . .

With Explicit Euler the solution approaches sin t as t ,

1
as the exact solution does, only if 0 < h <
25
This means h must be kept small, not to keep errors small,
but to have a stable numerical solution

The trapezoidal rule solution tends to sin t as t , as the

exact solution does, for every h > 0

This means h must be kept small only to keep errors small

Numerical Methods for Differential Equations p. 50/52

 Stiffness. . .

With Explicit Euler the solution approaches sin t as t ,

1
as the exact solution does, only if 0 < h <
25
This means h must be kept small, not to keep errors small,
but to have a stable numerical solution

The trapezoidal rule solution tends to sin t as t , as the

exact solution does, for every h > 0

This means h must be kept small only to keep errors small

For A-stable methods, there is no stability restriction on h

The step size is only restricted by accuracy

Numerical Methods for Differential Equations p. 50/52

 Stiffness. . .

Stiff differential equations are characterized by

homogeneous solutions being strongly damped

Example y = (y sin t) + cos t with 1

Stability regions of explicit methods are bounded, and

h D puts a strong stability restriction on h

Example Explicit Euler applied to the problem above

requires h 2/ 1

Numerical Methods for Differential Equations p. 51/52

 Stiffness. . .

Implicit methods with unbounded stability regions put no

restrictions on h, and the stepsize is only restricted by
accuracy requirement!

Example The implicit Euler applied to the problem above

requires only that

h2 y h2
ln sin tn
2 2
is sufficiently small, independently of

Numerical Methods for Differential Equations p. 52/52