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
.