7Implementation of ODE methods

IB Numerical Analysis

7.2 Solving for implicit methods

Implicit methods are often more likely to preserve nice properties like conservation

of energy. Since we have more computational power nowadays, it is often

preferable to use these more complicated methods. When using these implicit

methods, we have to come up with some way to solve the equations involved.

As an example, we consider the backward Euler method

y

n+1

= y

n

+ hf(t

n+1

, y

n+1

).

There are two ways to solve this equation for

y

n+1

. The simplest method is

functional iteration. As the name suggests, this method is iterative. So we use

superscripts to denote the iterates. In this case, we use the formula

y

(k+1)

n+1

= y

n

+ hf(t

n+1

, y

k

n+1

).

To do this, we need to find a

y

(0)

n+1

. Usually, we start

y

(0)

n+1

=

y

n

. Even better,

we can use some simpler explicit method to obtain our first guess of y

(0)

n+1

.

The question, of course, is whether this converges. Fortunately, this converges

to a locally unique solution if

λh

is sufficiently small, where

λ

is the Lipschitz

constant of

f

. For the backward Euler, we will require

λh <

1. This relies on

the contraction mapping theorem, which you may have met in IB Analysis II.

Does this matter? Sometimes it does. Usually, we pick an h using accuracy

considerations, picking the largest possible

h

that still gives us the desired

accuracy. However, if we use this method, we might need to pick a much smaller

h

in order for this to work. This will require us to compute much more iterations,

and can take a lot of time.

An alternative is Newton’s method. This is given by the formula

(I − hJ

(k)

)z

(k)

= y

(k)

n+1

− (y

n

+ hf(t

n+1

, y

(k)

n+1

))

y

(k+1)

n+1

= y

(k)

n

− z

(k)

,

where J

(k)

is the Jacobian matrix

J

(k)

= ∇f(t

n+1

, y

(k)

n+1

) ∈ R

N×N

.

This requires us to solve for

z

in the first equation, but this is a linear system,

which we have some efficient methods for solving.

There are several variants to Newton’s method. This is the full Newton’s

method, where we re-compute the Jacobian in every iteration. It is also possible

to just use the same Jacobian over and over again. There are some speed gains

in solving the equation, but then we will need more iterations before we can get

our y

n+1

.