INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS

I General method to construct a numerical integrator of ODEs 1

II The Forward Euler method 4

III The Backward Euler method 8

IV The Trapezoidal method 12

V Runge-Kutta methods 14

VI Higher-order Ordinary Differential Equations 16

VII The Verlet algorithm for second-order Ordinary Differential Equations 17

VIII Boundary Value Problems 21

In this chapter, we discuss the integration of Ordinary Differential Equations (ODEs).

Color code for the chapter:

▶ Definitions and important results are given in pink boxes.

▶ Skeletons of algorithms are in purple boxes.

▶ Examples are given in green boxes.

▶ Remarks and notions that can be skipped for the first read are in black boxes.

▶ Links are emphasized in blue.

Before starting, we want to stress that there is a full zoo of methods to integrate ODEs (some of them are
implemented in Python, see the documentation page of [Link]). It is impossible and useless to list
them all. We will present some of them which are particularly relevant to understand key concepts (convergence
of an integrator, consistency of an integrator, symplectic integrator, explicit or implicit integrator, stability of an
integrator, etc.) and which are heavily used in Physics.

We start by discussing the case of first-order ODEs and Initial Value Problems (IVPs), namely,
 #»
dx

 = f ( #»
x , t),
dt (1)
 #»
x (t ) = x#»,

i i

where #»x ∈ Rd is a d-dimensional vector, and where f : Rd × R → Rd . To be clear, #» x = (x1 , x2 , . . . , xd )

and f ( #»
x , t) = (f1 ( #»
x , t), f2 ( #»
x , t) . . . , fd ( #»
x , t)). The case of ODEs of larger order, as well as Boundary Value
Problems (BVPs) will be discussed at the end of this chapter.

1
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

I. General method to construct a numerical integrator of ODEs

1. Introduction

The basic idea behind the integration of ODEs is to discretize in time the Cauchy problem (1). Imagine that
you want to integrate in the time interval [ti , tf ] with N steps of size h = (tf − ti )/(N − 1). Your objective
is to compute x#»1 (h) , . . . , x# »
n
# »(h) such that x# »(h) is a good approximation of the exact solution x#»(t)
(h) , . . . , x
N n e
#»
evaluated at time tn = ti +nh, xe (tn ). Different algorithms come from different choices of the time discretization
of the derivative. This is discussed in the next paragraph.

2. Discretization of the time derivative

a. One-step methods

The idea is to perform Taylor expansion of the exact solution, and then to assume that the numerical integrator
is convergent, namely, that it approximates well the solution (see below for a proper definition). For instance, a
first-order Taylor expansion gives
dx#»e 1 d2 x#»e
x#»e (tn + h) = x#»e (tn ) + h (tn ) + h2 2 (tn ) + O(h3 ) = x#»e (tn ) + hf (x#»e (tn ), tn ) + O(h2 ). (2)
dt 2 dt
Then, if we assume that the method is convergent, x#»(t + h) = x# »(h) and x#»(t ) = x# »(h) and we eventually
e n n+1 e n n
get the following time discretization of Eq. (1):
x# n+1
»(h) = x# »(h) + hf (x# », t ).
n n n

This is called the Forward Euler method because the function f (i.e., the slope of the solution) is evaluated
at the previous time step (see Fig. 1). This method will be discussed below.

But we could have done the Taylor expansion the other way around, namely,
dx#»e
x#»e (tn ) = x#»e (tn + h) − h (tn + h) + O(h2 ) = x#»e (tn + h) − hf (x#»e (tn + h), tn + h) + O(h2 ), (3)
dt
leading to the time discretization
x# n+1
»(h) = x# »(h) + hf (x# », t + h).
n n+1 n

This is called the Backward Euler method because the function f is evaluated at the next time step (see
Fig. 1). This method wil also be discussed below.

We could also consider a point between tn and tn + h as a reference to perform the Taylor expansion, for instance
the midpoint tn + h/2. Then, the Taylor expansions are
h dx#»e 1 h 2 d2 x#»e
 
x#»e (tn ) = x#»e (tn + h/2) − (tn + h/2) + (tn + h/2) + O(h3 )
2 dt 2 2 dt2
and
h dx#»e 1 h 2 d2 x#»e
 
#» #»
xe (tn + h) = xe (tn + h/2) + (tn + h/2) + (tn + h/2) + O(h3 ),
2 dt 2 2 dt2
such that
dx#»e
x#»e (tn + h) = x#»e (tn ) + h (tn + h/2) + O(h3 ). (4)
dt
However, the value of the derivative at time tn + h/2 is not known so we need to express it as a function of the
values at times tn and tn + h, thanks to two other Taylor expansions:
 #»
dxe dx#»e h d2 x#»e

 (t n ) = (tn + h/2) − (tn + h/2) + O(h2 ),
dt dt 2 dt2

 dx#»e dx#»e h d2 x#»e
(tn + h/2) + O(h2 ),

 (tn + h) = (tn + h/2) +
dt dt 2 dt2

2
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

x
h/2 h/2

Forward Euler method

Exact Solution
Backward Euler method

(h)
Trapezoidal method
xn
t
tn tn+1

Figure 1: Illustration of three different methods to integrate an ODE. We show the behavior for one time
step of the Forward Euler method (in green), the Backward Euler method (in pink) and the Trapezoidal method
(in purple). The exact solution is represented as a black continuous line, and the dashed lines mark the tangent
lines to the exact solution at times tn and tn+1 . The black point marks the numerical solution at step n, and
the different colored points mark the numerical solution at step n + 1 for the different integration methods.

such that
dx#»e 1 dx#»e dx#»e
 
(tn + h/2) = (tn ) + (tn + h) + O(h2 ). (5)
dt 2 dt dt
We eventually obtain the following time discretization:
»(h) = x# »(h) + h f (x# »(h) , t ) + f (x# »(h) , t + h) .
x# n+1
h i
n n n n+1 n
2
This is called the Trapezoidal method because the function f is evaluated at the previous time step and at the
next time step (see Fig. 1). This method will be discussed below.

b. Multi-step methods

So far, we have considered time discretizations which only involve the knowledge of the solution at the previous
step. These methods are called one-step methods. However, there are schemes which require the knowledge of
the solution at the two previous steps or more (multi-step methods). We provide below an example. We again
start from the Taylor expansion of x#»e at time tn + h, see Eq. (2). We then use a Taylor expansion of the time
derivative of x#»e at time tn − h to express the second derivative:
dx#»e dx#»e d2 x#»e
(tn − h) = (tn ) − h 2 (tn ) + O(h2 )
dt dt dt
d2 x#»e
=⇒ h 2 (tn ) = f (x#»e (tn ), tn ) − f (x#»e (tn − h), tn − h) + O(h2 ).
dt
This eventually leads to
3h #» h
x#»e (tn + h) = x#»e (tn ) + f (xe (tn ), tn ) − f (x#»e (tn − h), tn − h) + O(h3 ),
2 2
from which we deduce the following time discretization:
»(h) = x# »(h) + h 3f (x# »(h) , t ) − f (x# »(h) , t − h) .
x# n+1
h i
n n n n−1 n
2
This is called the Adams-Bashforth method.

c. Conclusion

Depending on the time discretization for the derivatives, you can construct a large family of numerical integrators
of ODEs. The starting point is always a Taylor expansion of the exact solution at different time steps (your
creativity is the limit!).

3
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

General definition of a numerical integrator

A numerical integrator of an ODE of the form given by Eq. (1) is a relation giving the numerical estimate
of the solution at step n + 1 as a function of the numerical estimates of the solution in the present and in
the k past steps:
»(h) , x# »(h) , x# »(h) , . . . , x# »(h) , t ) = #»
G(x# n+1 0. (6)
n n−1 n−k n

The method is explicit if the numerical estimate of the solution at step n + 1 only depends on the numerical
estimates in the past. It can then be written as

x# n+1
» = H(x# »(h) , x# »(h) , . . . , x# »(h) , t )
n n−1 n−k n (7)

for a k-step method.

The method is implicit if the numerical estimate of the solution at step n depends on the numerical estimates
in the past and in the present.

II. The Forward Euler method

In this section, we discuss the Forward Euler method introduced above.

1. Definition

The Forward Euler method is defined as follows.

Definition of the Forward Euler method

For the IVP given by Eq. (1), the Forward Euler method is a numerical integrator of ODEs given by the
following recurrence relation:
x# n+1
»(h) = x# » + hf (x# », t ),
n n n (8)
for tn = ti + nh and h the step size.

The Forward Euler method is explicit. Therefore, we can propose a simple algorithmic representation.

Algorithm for the Forward Euler method

x ← xi ▷ Initial condition.
t ← ti ▷ Initial time.
store(x) ▷ Store the initial value of x.
store(t) ▷ Store the initial value of t.
while t < tf do
x ← x + hf (x, t) ▷ Update the solution.
t←t+h ▷ Update the time.
store(x) ▷ Store the current value of x.
store(t) ▷ Store the current value of t.

We discuss the accuracy of the method in the next section.

2. Local truncation error and consistency

We come back to the Taylor series at the root of the Forward Euler method, see Eq. (2). You can see that we
have truncated all terms of order h2 and above to define the Forward Euler method, see Eq. (8). We thus say
that the local truncation error is O(h2 ).

Local truncation error

The local truncation error τn at step n is defined as the difference between x# » n
(h) (estimated numerical

value of the solution at time tn ) and x#»e (tn ) (exact value of the solution at time tn ), assuming that no error

4
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

was made in the previous steps. In other words,

τn = ∥x# »
n
(h)
− x#»e (tn )∥ assuming x# »k (h) = x#»e (tk ) ∀k < n. (9)

From the local truncation error, we define the consistency of a numerical integrator.

Consistency of a numerical integrator

A numerical integrator of ODEs is consistent if the time discretization approximates the ODE well. Said
differently, the exact solution to an ODE should verify the recurrence relation up to terms of order hq with
h the step size and q > 1 in the limit h → 0.

Mathematically, a numerical integration is consistent if its local truncation error τn goes to 0 when the step
size h goes to 0 at least as h1 , or equivalently, if
τn
∀n ∈ J1, N K lim = 0. (10)
h→0 h

From the above definitions, we can deduce some properties of the Forward Euler method.

Consistency of the Forward Euler method

The Forward Euler method has a local truncation error τn = O(h2 ): it is therefore consistent.

We end this section by stressing that the consistency does not tell you whether the numerical integrator
gives you a good estimate of the solution of an ODE. Instead, it just tells you that the Cauchy problem is
well captured by the numerical integrator. Whether the numerical integrator gives a reasonable estimate of the
solution is called convergence, and is discussed in the next section.

3. Global truncation error and convergence

Of course, when you integrate an ODE on a interval [ti , tf ], local truncation errors accumulate. Therefore, a
better assessment of the accuracy of a numerical integrator is given by looking at the global truncation error.

Global truncation error

The global truncation error en at step n is defined as the difference between x# » n
(h) (estimated numerical

value of the solution at time tn ) and x#»e (tn ) (exact value of the solution at time tn ), taking into account all
possible errors in the previous steps. In other words,

en = ∥x# »
n
(h)
− x#»e (tn )∥. (11)

From the global truncation error, we define the convergence of a numerical integrator.

Convergence of a numerical integrator

A numerical integrator of ODEs is convergent if the numerical solution approaches the exact solution at
any time when the step size h goes to 0.

Mathematically, a numerical integrator is convergent if

∀t ∈ [ti , tf ] lim en = 0. (12)

h→0
n→+∞
nh=t−ti

The way en goes to 0 when we take the two limits h → 0 and n → +∞ but with nh < +∞ makes it possible
to define the order of the numerical integrator: the method is of order p (with p > 0) if en = O(hp ) when
h → 0, n → +∞, nh < +∞.

The convergence property is, of course, a property that you absolutely need in order to integrate an ODE
in Physics.

5
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

We now illustrate the convergence property on the Forward Euler method. It is possible to prove that the Forward
Euler method is of order 1, see D. F. Griffiths and D. J. Higham, Numerical Methods for Ordinary Differential
Equations. However, the proof is complicated, and we just give an example below.
Imagine that you want to solve numerically the Cauchy problem
dx
= x, x(0) = 1,
dt
on [0, 1]. We know that the exact solution is xe (t) = et .
(h) (h) (h) (h)
Now, let’s look at the Forward Euler method: xn+1 = xn + hxn = xn (1 + h). This is a geometric progression,
(h)
and this implies that xn = (1 + h)n . From this, we can compute the global truncation error:

en = (1 + h)n − enh = en ln(1+h) − enh .

We now perform the two limits h → 0, n → +∞ but with nh = t finite and we obtain
 
2
/2+O(h3 )] 2
/2+nO(h3 ) 2 th
en = en[h−h − enh = enh e−nh − 1 = et e−th/2+tO(h ) − 1 = et + tO(h2 )
2
1 t
= te h + higher-order terms.
2
This implies that en = O(h) as mentioned above.

From the above definitions, we can deduce some properties of the Forward Euler method.
Convergence of the Forward Euler method
The Forward Euler method has a global error O(h): it is of order 1 and it is therefore convergent.

Being of order 1, to increase by one digit the accuracy of the Forward Euler method (i.e., to divide by a
factor of 10 the global truncation error), you have to divide h by a factor of 10: you have to make 10 times
more steps and your computation time is also multiplied by a factor of 10.

For the majority of numerical integrators of ODEs (see though a counter-example below), if the local truncation error
is O(hp+1 ), the global error is O(hp ). Therefore, looking at the local truncation error is most of the time enough to
decide whether the method is convergent.

4. Stability and stiffness

In the previous section, we have focused on the convergence of a numerical integrator, namely, with what happens
when h → 0. In practice, to reduce the computation time, you want to increase the value of h. Therefore, a
natural question which emerges is how large you can make the step size to get a result that is still a correct
estimate of the exact solution. For most methods, you cannot take h too large because numerical integrators
of ODEs can become unstable (the result of the integration diverges from the exact solution). The theory of
the stability of numerical integrators brings answers to the above question. However, it is relatively complex.
Therefore, we will not delve into the details of such a theory, but instead illustrate stability issues when numerically
integrating an ODE in the example below.
Consider the Cauchy problem
dx
= −ax, x(0) = 1,
dt
on [0, 1] with a > 0. The exact solution is xe (t) = e−at , while the Forward Euler method leads to
(h)
xn+1 = x(h) (h) (h)
n − ahxn = xn (1 − ah).

(h)
This is a geometric progression and for any n we have xn = (1 − ah)n . The exact solution of the ODE xe (t) = e−at

6
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

(h)
decays exponentially, and xn should be a decaying sequence with n. However, this is only the case if −1 < 1−ah < 1,
otherwise the sequence diverges. The two behaviors are shown Fig. 2. Therefore, the integrator is stable if h < 2/a.

For certain ODEs, which are called stiff, you need to consider very small values of h in order for the method to
be stable.

Stiff ODEs
Some ODEs are stiff, meaning that they may require a very small time step h in order for the numerical
integrator to be stable. This often occurs when you have multiple characteristic timescales with different
orders of magnitude in the differential equation.

The value of the step size h to recover stability depends on the ODE and also on the numerical solver.

h < 2/a h > 2/a

1.2 3
e−at
1 (h) 2
xn

0.8 1

0.6 0
x

0.4 −1

0.2 −2 e−at
(h)
xn
0 −3
0 2 4 6 8 10 0 2 4 6 8 10
t t

Figure 2: Stability issues of the Forward Euler method. We show the exact solution of the ODE dx/dt = −ax
(h)
for a > 0 and x(0) = 1 along with the numerical solution xn obtained from the Forward Euler method. In the
left panel, the step size is h < 2/a so that the integration scheme is stable. In the right panel, the step size is
h > 2/a and the integration scheme is unstable: the numerical solution diverges.

5. Conclusion

We can summarize the advantages and drawbacks of the Forward Euler method.

Advantages and Drawbacks of the Forward Euler method

Advantages:

▶ The Forward Euler method is fast (few operations at each time step).

▶ The Forward Euler method is simple to implement and only requires to store the value at the previous
step.

Drawbacks:

▶ The Forward Euler method is not very precise.

▶ The Forward Euler method can be unstable.

7
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

III. The Backward Euler method

In this section, we discuss the Backward Euler method introduced above.

1. Definition

The Backward Euler method is defined as follows

Definition of the Backward Euler method

For the IVP given by Eq. (1), the Backward Euler method is a numerical integrator of ODEs given by the
following recurrence relation:
x# n+1
»(h) = x# » + hf (x# », t
n n+1 n+1 ), (13)
for tn = ti + nh and h the step size.

The Backward Euler method is a priori implicit because x# n+1

»(h) also appears in the right-hand side of Eq. (13).
There are three strategies to implement the algorithm.

First strategy: making the Backward Euler method explicit. In some cases, you can manipulate Eq. (13)
»(h) = g(x# »(h) , t
analytically in order to obtain an explicit expression of x# n+1 n n+1 ), see the example below. Once
you have found an explicit expression, you can implement it in a similar way as the Forward Euler method.

Algorithm for the Backward Euler method when made explicit

x ← xi ▷ Initial condition.
t ← ti ▷ Initial time.
store(x) ▷ Store the initial value of x.
store(t) ▷ Store the initial value of t.
while t < tf do
t←t+h ▷ Update the time.
x ← x + hg(x, t) ▷ Update the solution.
store(x) ▷ Store the current value of x.
store(t) ▷ Store the current value of t.

We present an example where we can turn the implicit equation defining the Backward Euler method into an explicit
expression. Consider the Cauchy problem already introduced in the example above
dx
= −ax, x(0) = 1,
dt
(h) (h) (h)
on [0, 1] with a > 0. The Backward Euler method is defined by the equation xn+1 = xn − ahxn+1 , that you can
(h)
transform into an explicit equation for xn+1 :
(h)
(h) xn
xn+1 = .
1 + ah

Of course, it is not always possible to make the Backward Euler method explicit. We thus propose two other
strategies.

Second strategy: finding the solution by iteration. Equation (13) implies that x# n+1 »(h) is a fixed point
#» # »(h) #»
of the function F : x 7→ xn + hf ( x , tn+1 ). There is a mathematical theorem which tells you that, under
reasonable properties of the function F , if you consider the sequence (y#»p )p∈N defined by the recursion relation
y# p+1
» = F (y#»), this sequence converges to the fixed point of F , which is precisely the solution x# »(h) of the
p n+1
Backward Euler equation. As a consequence, the idea is to iterate the function F , starting from y#»0 = x# »
n
(h) until

the sequence (y#»p )p∈N converges. In practice, you have to stop the iteration when the relative change between
two iterations is smaller than a given threshold ϵ.

8
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

Algorithm for the Backward Euler method with an iteration scheme

x ← xi ▷ Initial condition.
t ← ti ▷ Initial time.
store(x) ▷ Store the initial value of x.
store(t) ▷ Store the initial value of t.
while t < tf do
t←t+h ▷ Update the time.
y←x ▷ Initiate the sequence for the research of the fixed point.
z ← +∞ ▷ z represents the value of y at the previous iteration: this value is stored to check whether
y converges to the fixed point.
while ∥y − z∥ > ϵ∥y∥ do ▷ Iterate the function F until the sequence converges to the fixed point up
to a given threshold ϵ.
z←y
y ← x + hf (y, t)
x←y ▷ Update the value of the solution.
store(x) ▷ Store the current value of x.
store(t) ▷ Store the current value of t.

However, this iterative method can be unstable: the success of the iterative process depends on the mathematical
properties of F , and therefore on the value of the step size h which should be small enough. We illustrate this
in the example below, and we provide next a last strategy to solve Eq. (13).

We come back to the Cauchy problem considered above

dx
= −ax, x(0) = 1,
dt
on [0, 1] with a > 0, and we apply the iteration scheme at step n + 1. We thus define the sequence (yp )p∈N such
(h) (h)
that y0 = xn and yp+1 = xn − ahyp . You can thus obtain a generic expression for yp (this is a good mathematical
exercise):
(h)
xn 
1 − (−ah)p+1 .

yp =
1 + ah
(h)
You then see that (yp )p∈N converges to the fixed point xn /(1 + ah) (see the previous example where we made the
implicit Euler method explicit) only if −1 < −ah < 1, i.e., h < 1/a. We illustrate the two behaviors h < 1/a and
h > 1/a in Fig. 3. You can note that this criterion is roughly similar to the stability criterion for the Forward Euler
method (we found h < 2/a in that case).

Third strategy: finding the solution with a root-finding algorithm. Equation (13) implies that x# n+1 »(h) is a
#» #» # » (h) #» # » (h) #»
root of the function G : x 7→ x − xn − hf ( x , tn+1 ), namely, G(xn+1 ) = 0 . There is a well-known method
to find the root of a function, which is called the Newton method, that we detail now. The idea is again to build
a sequence (y#»p )p∈N with y#»0 = x# »
n
#»)
(h) , and to iterate such that (y
p p∈N converges to the root of G. The iteration
# » #» #» #»
is defined as follows: yp+1 = yp − δp where δp = (δp,1 , . . . , δp,d ) ∈ Rd is the solution of the linear system

∂G1 #» ∂G1 #»
 
(y ) . . . (yp )
 ∂y1 p ∂yd  #»

.. ..  δp = G(y#»p ),


 . .  (14)
 ∂Gd ∂Gd #» 
(y#»p ) . . . (yp )
∂y1 ∂yd

with G( #»
y ) = (G1 ( #»
y ), . . . , Gd ( #»
y )). The matrix of the partial derivatives in the left hand-side is called the
Jacobian matrix. This method is implemented in Python, look at the documentation page of [Link].
As for the iteration procedure, you have to truncate the sequence (y#»p )p∈N when the relative change between
two iterations is smaller than a given threshold ϵ. We illustrate the method with an example below, provide a
mathematical proof in a remark, and finally give a schematic algorithm.

9
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

h < 1/a h > 1/a

4
2 (h)
F (y) = xn − ahy
y
2
1.5
F (y), y

F (y), y
0
1

0.5 −2
(h)
F (y) = xn − ahy
y
0 −4
0 0.5 1 1.5 2 −2 −1 0 1 2
y y

Figure 3: Convergence issue when solving the Backward Euler equation with an iteration scheme. We
illustrate the method to solve for step n + 1 the ODE dx/dt = −ax for a > 0 from step n using the Backward
(h)
Euler method with an iteration scheme. We define the sequence y0 = xn (violet disks) and yp+1 = F (yp ) with
(h)
F (y) = xn − ahy. For a step size h small enough (h < 1/a), the iteration scheme converges to the fixed point
(violet square). Instead for a too large step size h (h > 1/a), the iteration scheme diverges and is unstable.

We come back to the Cauchy problem considered above

dx
= −ax, x(0) = 1,
dt
on [0, 1] with a > 0, and we apply the root-finding scheme at step n + 1. We thus define the sequence (yp )p∈N such
(h)
that y0 = xn and yp+1 = yp − δp with

G(yp )
δp = , with G(y) = y − xn(h) + ahy.
G′ (yp )
(h) (h)
You then easily get that δp = yp − xn /(1 + ah) so that yp+1 = xn /(1 + ah). As a result, the sequence always
(h)
converges to xn+1 (see the above example where we made the implicit Backward Euler method explicit) in one step,
whatever the value of h.

In this remark, we prove the Newton method. Imagine that y#»p is close to x# n+1 »(h) . The difference δ#» = y#» − x# »(h)
p p n+1
#»
is thus much smaller than y#»p : ∥δp ∥ ≪ ∥y#»p ∥. If we use the fact that x# n+1
»(h) is a root of G, it means that

∀ i ∈ J1, dK Gi (x# n+1

»(h) ) = 0.

»(h) = y#» − δ#», we have that

If we now use the fact that x# n+1 p p

∀ i ∈ J1, dK Gi (yp,1 − δp,1 , . . . , yp,d − δp,d ) = 0.

We can then do a Taylor expansion at first order, and we obtain that

d
∂Gi
∀ i ∈ J1, dK Gi (y#»p ) − (y#»p )δp,j = 0.
X

j=1
∂yj

Putting this into a matrix form gives you Eq. (14).

10
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

Algorithm for the Backward Euler method with a root-finding scheme

x ← xi ▷ Initial condition.
t ← ti ▷ Initial time.
store(x) ▷ Store the initial value of x.
store(t) ▷ Store the initial value of t.
while t < tf do
t←t+h ▷ Update the time.
y←x ▷ Initiate the sequence for the research of the root.
z ← +∞ ▷ z represents the value of y at the previous iteration: this value is stored to check whether
y converges to the root.
while ∥y − z∥ > ϵ∥y∥ do ▷ Iterate until convergence to the root up to a given threshold ϵ.
z←y
Geval ← G(y) ▷ Evaluate G(y#»p ).
Jeval ← J(y) #»
▷ Evaluate the Jacobian matrix J(yp ) of the partial derivatives.
#»
δ ← solve_linear_system(Jeval , Geval ) ▷ Solve the linear system defining δp .
y ←y−δ ▷ Update the value of the sequence.
x←y ▷ Update the value of the solution.
store(x) ▷ Store the current value of x.
store(t) ▷ Store the current value of t.

We discuss the properties of the Backward Euler method in the next section.

2. Properties

We give the main properties of the Backward Euler method, which are the consequences of the Taylor expansion
given by Eq. (3).

Consistency and convergence of the Backward Euler method

The Backward Euler method has a local truncation error O(h2 ): it is consistent.

The Backward Euler method has a global truncation error O(h) (method of order 1): it is convergent.

Its characteristics are thus similar to the Forward Euler method. However, it may require much more operations
than the Forward Euler method if you use the strategy of the fixed point or the one of root-finding to solve
Eq. (13). Still, the Backward Euler method has one advantage with respect to the Forward Euler method: it is
more stable if you can make the Backward Euler method explicit, or if you use a root-finding algorithm. We
illustrate this in the example below.

We come back to the Cauchy problem

dx
= −ax, x(0) = 1,
dt
on [0, 1] with a > 0. We have seen above that the Backward Euler method can be made explicit:
(h)
(h) xn
xn+1 = .
1 + ah
(h)
This implies that xn has a geometric progression, leading to
 n
1 1
x(h)
n = = .
(1 + ah)n 1 + ah

For all values of h, this sequence decays with n, because 0 < 1/(1 + ah) < 1, and the method is thus stable.

We have also discussed above the stability conditions of the iteration scheme and of the root-finding scheme.

We end this section by summarizing the advantages and drawbacks of the Backward Euler method.

11
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

Advantages and Drawbacks of the Backward Euler method

Advantages:

▶ The Backward Euler method is relatively simple to implement and only requires to store the value
at the previous step.

▶ The Backward Euler method is more stable if you can make it explicit or if you use a root-finding
algorithm: it is better for stiff problems.

Drawbacks:

▶ The Backward Euler method may be slow if you have to iterate in order to solve the equation defining
the value of the solution at the next step.

▶ The Backward Euler method is not very precise.

IV. The Trapezoidal method

In this section, we briefly discuss the Trapezoidal method introduced above.

1. Definition

The Trapezoidal method is defined as follows

Definition of the Trapezoidal method

For the IVP given by Eq. (1), the Trapezoidal method is a numerical integrator of ODEs given by the
following recurrence relation:

»(h) = x# » + h [f (x# », t ) + f (x# », t

x# n+1 (15)
n n n n+1 n+1 )] ,
2
for tn = ti + nh and h the step size.

(h)
The Trapezoidal method is a priori implicit because xn+1 also appears in the right-hand side of Eq. (15). There
are four strategies to implement the algorithm.

First strategy: making the Trapezoidal method explicit. As for the Backward Euler method, you can
manipulate Eq. (15) analytically in order to obtain an explicit expression of x# n+1
»(h) = g(x# »(h) , t , t
n n n+1 ), see the
example below. Once you have found an explicit expression, you can implement it in a similar way as the Forward
Euler method.
We present an example where we can turn the implicit equation defining the Trapezoidal method into an explicit
expression. Consider again the Cauchy problem
dx
= −ax, x(0) = 1,
dt
on [0, 1] with a > 0. The Trapezoidal method is defined by the equation

(h) ah h (h) (h)

i
xn+1 = x(h)
n − xn + xn+1 ,
2
(h)
that you can transform into an explicit equation for xn+1 :

(h) 1 − ah/2
xn+1 = x(h)
n .
1 + ah/2

Of course, as before, this is not always possible to make the Trapezoidal method explicit. In the latter case,
you have to use the methods which have been introduced in the previous section about the Backward Euler

12
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

method (second strategy: finding the solution by iteration and third strategy: finding the solution with
a root-finding algorithm). Or you can use a much simpler method which is specific to the Trapezoidal method.

Fourth strategy: using a predictor-corrector strategy. The idea is to make a first guess of x# n+1»(h) (predictor
#»
y ) using a Forward Euler method, and then to use this predictor value to compute the actual value of x# n+1»(h)
#
(corrector) by replacing xn+1 » (h) #»
by y in the right-hand side of Eq. (15). In other words, you have turned an
implicit one-stage method into a two-stage explicit method. The combination of the Trapezoidal method with a
predictor-corrector strategy is sometimes called the Heun’s method.

Algorithm for the Trapezoidal method with the predictor-corrector scheme

x ← xi ▷ Initial condition.
t ← ti ▷ Initial time.
store(x) ▷ Store the initial value of x.
store(t) ▷ Store the initial value of t.
while t < tf do
t′ ← t ▷ Store the time at step n.
t←t+h ▷ Update the time at step n + 1.
y ← x + hf (x, t′ ) ▷ Compute the predictor.
x ← x + (h/2)f (x, t′ ) + (h/2)f (y, t) ▷ Compute the corrector.
store(x) ▷ Store the current value of x.
store(t) ▷ Store the current value of t.

Another way of seeing the predictor-corrector strategy is to notice that it is equivalent to the iteration scheme
(second strategy) but truncated after the first iteration. As a consequence, this method can also be unstable for
stiff problems, as illustrated in the example below.

We come back to the Cauchy problem considered above

dx
= −ax, x(0) = 1,
dt
(h) (h)
on [0, 1] with a > 0, and we apply the predictor-corrector scheme at step n + 1. We thus have y = xn − ahxn
and then
(h) ah (h) ah
xn+1 = xn(h) − x − y,
2 n 2
leading to
(ah)2
 
(h)
xn+1 = x(h)
n 1 − ah + .
2
(h)
You can recognize the Taylor expansion of t 7→ e−at (the exact solution) up to second order. The estimates xn
(h)
follow a geometric progression and we can eventually obtain a closed expression for xn :
n
(ah)2

x(h)
n = 1 − ah + .
2
(h)
The exact solution is an exponentially decaying function, and xn should be a decaying sequence with n. However,
this is only the case if −1 < 1 − ah + (ah)2 /2 < 1, resulting in h < 2/a in order for the method to be stable.
Otherwise, the sequence diverges. This is the same stability criterion as for the Forward Euler method.

We discuss the properties of the Trapezoidal method in the next section.

2. Properties

We give the main properties of the Trapezoidal method, which are the consequences of the Taylor expansions
given by Eq. (4) and Eq. (5).

13
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

Consistency and convergence of the Trapezoidal method

The Trapezoidal method has a local truncation error O(h3 ): it is consistent.

The Trapezoidal method has a global truncation error O(h2 ) (method of order 2): it is convergent. Said
differently, to increase by one digit the accuracy of the
√method (i.e., to divide by a factor of 10 the global
truncation error), you have to divide h by a factor of 10 ≃ 3 only: you have to make 3 times more steps
and your computation time is also multiplied by a factor of 3.
We end this section by summarizing the advantages and drawbacks of the Trapezoidal method.
Advantages and Drawbacks of the Trapezoidal method
Advantages:

▶ The Trapezoidal method is relatively simple to implement and only requires to store the value at the
previous step.

▶ The Trapezoidal method is more stable if you can make it explicit or if you use the Newton method:
it is better for stiff problems.

▶ The Trapezoidal method is more precise than the Euler methods (see also Fig. 1).

Drawbacks:

▶ The Trapezoidal method may be slow if you have to iterate in order to solve the equation defining
the value of the solution at the next step.

V. Runge-Kutta methods
We now discuss Runge-Kutta methods which are heavily used in Physics for their high degree of precision and
their easy implementation.

1. Definition

In the previous sections, we have discussed the Euler methods and the Trapezoidal method. The latter is more
precise than the formers. The change in precision comes from a different choice of points from which the Taylor
expansions in Eqs. (2), (3), (4) and (5) are performed: the Taylor expansion is done from t = tn for the Forward
Euler method, t = tn+1 for the Backward Euler method, and t = tn + h/2 for the Trapezoidal method. To derive
Runge-Kutta methods, you just have to consider several points between tn and tn+1 and play with the Taylor
expansions to cancel as many powers of h as you want in order to increase the accuracy. You can construct an
entire family of Runge-Kutta methods, see J. C. Butcher, Numerical Methods for Ordinary Differential Equations.
In these notes, we only discuss the most famous Runge-Kutta method: RK4.
Definition of the RK4 method
For the IVP given by Eq. (1), the RK4 method is a numerical integrator of ODEs given by the following
recurrence relations:  #»
k1 = hf (x# »n , tn ),
#» #»

# »


k2 = hf (xn + k1 /2, tn + h/2),


 #» #»

k3 = hf (x# »

n + k2 /2, tn + h/2), (16)
#» # » + k#», t + h),
k = hf ( x


 4 n 3 n
»(h) = x# »(h) + 1 k#» + 2k#» + 2k#» + k#» ,

x# n+1

  

n 1 2 3 4
6
for tn = ti + nh and h the step size.

The RK4 method is explicit. It is a one-step method (because you only require the knowledge of the solution at
the previous step) but it has multiple stages. We can propose a simple algorithmic representation. Note however
that this method is already implemented in Python, see the documentation.

14
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

Algorithm for the RK4 method

x ← xi ▷ Initial condition.
t ← ti ▷ Initial time.
store(x) ▷ Store the initial value of x.
store(t) ▷ Store the initial value of t.
while t < tf do
k1 ← hf (x, t)
k2 ← hf (x + k1 /2, t + h/2)
k3 ← hf (x + k2 /2, t + h/2)
k4 ← hf (x + k3 , t + h)
x ← x + (k1 + 2k2 + 2k3 + k4 )/6 ▷ Update the solution.
t←t+h ▷ Update the time.
store(x) ▷ Store the current value of x.
store(t) ▷ Store the current value of t.

We discuss the properties of the RK4 method in the next section.

2. Properties

We give the main properties of the RK4 method, which are again the consequences of the Taylor expansions
which define it.
Consistency and convergence of the RK4 method
The RK4 method has a local truncation error O(h5 ): it is consistent.

The RK4 method has a global truncation error O(h4 ) (method of order 4): it is convergent. Said differently,
to increase by one digit the accuracy of the√method (i.e., to divide by a factor of 10 the global truncation
error), you have to divide h by a factor of 4 10 ≃ 1.8 only: you have to make less than twice more steps
and your computation time is also multiplied by a factor of 1.8 roughly.

Although the RK4 method is very precise, it can be unstable when trying to solve a stiff problem, as exemplified
below.
Consider the Cauchy problem
dx
= −ax, x(0) = 1,
dt
on [0, 1] with a > 0. The RK4 method applied to this problem gives:
 (h)

 k1 = −ahxn ,


  
ah

 (h)
k2 = −ah 1 − 2 xn ,




  
ah ah (h)
k3 = −ah 1 − 1− xn ,





 2 2

    

 ah ah (h)
k4 = −ah 1 − ah 1 −
 1− xn ,
2 2

such that you eventually get that

(ah)2 (ah)3 (ah)4

 
(h)
xn+1 = x(h)
n 1 − ah + − + .
2 6 24

In the brackets, you recognize the Taylor expansion of t 7→ e−at (which is the exact solution) up to fourth order.
This is a geometric progression and for any n we have
n
(ah)2 (ah)3 (ah)4

(h)
xn = 1 − ah + − + .
2 6 24

15
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

(h)
The exact solution is an exponentially decaying function, and xn should be a decaying sequence with n. However,
this is only the case if the common ratio (the quantity into brackets) lies between −1 and 1. You can prove that it
amounts to " 1/3 #
√ 1/3

1 2 2/3

ah < 4 − 10 √ +2 43 + 9 29 ≃ 2.78529...
3 43 + 9 29
Therefore, the RK4 method is stable if the above inequality is satisfied. Otherwise, the sequence diverges and the
RK4 method is unstable. You can note that the stability criterion derived above is less restrictive than for the Forward
Euler method (it was ah < 2 in that case).

We end this section by summarizing the advantages and drawbacks of the RK4 method.

Advantages and Drawbacks of the RK4 method

Advantages:

▶ The RK4 method is relatively simple to implement and only requires to store the value at the previous
step.

▶ The RK4 method is relatively fast (few operations at each time step).

▶ The RK4 method is much more precise than the Euler methods or the Trapezoidal method.

Drawbacks:

▶ The RK4 method can be unstable when you have to solve stiff problems. However, it is more stable
than the Forward Euler method.

To cure the stability issues of the RK4 methods, implicit Runge-Kutta methods have been introduced, but they are
not discussed in these lecture notes.

VI. Higher-order Ordinary Differential Equations

So far, we have only discussed Initial Value Problems (IVPs) with first-order differential equations. We now
discuss how to deal with higher-order differential equations in general. In the next section, we consider particular
second-order differential equations for which specific methods have been introduced.

Any ODE of order q > 1 can be recast into a system of q first-order ODEs by introducing auxiliary variable.
Imagine that you want to solve a scalar ODE of the form

dq x dq−1 x
 
dx
= f x, , . . . , q−1 , t ,
dtq dt dt

where x ∈ R. Then, you can define x0 = x, x1 = dx/dt, . . . , xq−1 = dq−1 x/dtq−1 such that the scalar ODE
of order q can be written as a system of q coupled first-order ODEs
dx0


 = x1 ,
dt






 dx1
= x2 ,


dt





..

. ,





 dxq−2

 = xq−1 ,
dt





 dxq−1 = f (x0 , x1 , . . . , xq−1 , t)



dt

16
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

that you can eventually recast in a vectorial form:

d #»
x
= F ( #»
x , t), with #»
x = (x0 , x1 , . . . , xq−1 ) and F ( #»
x , t) = (x1 , x2 , . . . , xq−1 , f (x0 , x1 , . . . , xq−1 , t)).
dt
You can then apply any of the methods introduced above to deal with first-order ODEs.

VII. The Verlet algorithm for second-order Ordinary Differential Equations

In this section, we discuss the numerical integration of Cauchy problems of the form:
 2 #»
d x


2
= f ( #»x ),
dt




#»
x (ti ) = x#»i ,

(17)

#»


 d x (t ) = v#»,



i i
dt
where #»
x ∈ Rd is a d-dimensional vector, and where f : Rd → Rd . Note that this Cauchy problem is not the
most general one for second-order differential equations for two reasons:

▶ the function in the right-hand side only depends on #»

x and not on its first derivative;

▶ the function in the right-hand side does not depend explicitely on time: the system is autonomous.

Despite these specificities, being able to solve such equations is of paramount importance because they naturally
emerge in Physics when solving Newton’s equations for the mechanics of conservative systems.
We can rewrite the above Cauchy problem as a system of vectorial first-order ODEs by introducing the auxiliary
variable #»
v = d #»
x /dt (the velocity):  #»
dx

 = #»
v,
dt



#»


 d v = f ( #»


x ),

dt (18)
#» #»




 x (ti ) = xi ,



 #»
v (t ) = v#»,


i i

1. Definition

We now define the velocity Verlet algorithm to solve the above Cauchy problem. For more details, you can
read M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids.

Definition of the velocity Verlet algorithm

For the IVP given by Eq. (17) [or Eq. (18)], the velocity Verlet algorithm is a numerical integrator given by
the following recurrence relations:
2

»(h) = x# »(h) + hv#»(h) + h f (x# »(h) ),
x# n+1

n n n
2

(19)
v# n+1
»(h) = v#»(h) + h
f (x# » # »(h) ) ,

(h) ) + f (x
  
n n n+1
2
for tn = ti + nh and h the step size.

The velocity Verlet method is explicit. We can thus propose a simple algorithmic representation.

17
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

Algorithm for the velocity Verlet method

x ← xi ▷ Initial position.
v ← vi ▷ Initial velocity.
a ← f (xi ) ▷ Initial acceleration.
t ← ti ▷ Initial time.
store(x) ▷ Store the initial value of position.
store(v) ▷ Store the initial value of velocity.
store(t) ▷ Store the initial value of time.
while t < tf do
x ← x + hv + (h2 /2)a ▷ Update the position.
v ← v + (h/2)a ▷ Update partially the velocity.
a ← f (x) ▷ Update the acceleration.
v ← v + (h/2)a ▷ Update the velocity.
t←t+h ▷ Update the time.
store(x) ▷ Store the current value of position.
store(v) ▷ Store the current value of velocity.
store(t) ▷ Store the current value of time.

We discuss the properties of the velocity Verlet method in the next section.

2. Properties

Consistency and convergence of the velocity Verlet method

The velocity Verlet method has a local truncation error O(h4 ): it is consistent.

The velocity Verlet method has a global truncation error O(h2 ) (method of order 2): it is convergent. Said
differently, to increase by one digit the accuracy of the
√ method (i.e., to divide by a factor of 10 the global
truncation error), you have to divide h by a factor of 10 ≃ 3 only.

The RK4 method is more precise than the velocity Verlet method. So why should we prefer the velocity Verlet
method when integrating Cauchy problems of the form given by Eq. (17)? We give the main motivation below
and illustrate it with an example after.

Symplecticness of the velocity Verlet algorithm

The velocity Verlet algorithm is symplectic. In particular, for conservative systems for which the total
energy is constant, the algorithm approximately preserves this conservation law.

Instead, the Runge-Kutta 4 method is not symplectic. For conservative systems, the total energy system-
atically decreases.

We consider the Cauchy problem corresponding to a harmonic oscillator:

 2
d x

 = −ω 2 x,
dt2





x(0) = 1,



 dx (0) = 0.



dt
We integrate it with the same step size h using the velocity Verlet algorithm and the Runge-Kutta 4 method. The
results are presented Fig. 4.

We first analyze the numerical solution obtained in a small time interval, and we confront it with the exact solution
cos(ωt) (see the left panels). Both methods agree well with the exact solution, although the numerical estimates
from the velocity Verlet algorithm start to deviate from the exact solution at the end of the time window. This comes

18
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

from the fact that the Runge-Kutta 4 method has a better global truncation error than the velocity Verlet algorithm.

We now analyze the numerical solution obtained in a larger time interval. The harmonic oscillator is a conservative
system, meaning that its total energy
 2
1 dx 1
E= m + mω 2 x2
2 dt 2
is constant during the dynamics. We now assess if this property is preserved by the velocity Verlet algorithm and
the Runge-Kutta 4 method. It turns out that for the velocity Verlet algorithm, the energy fluctuates around a value
which is close to the exact value of the energy. In other words, energy conservation is approximately preserved by the
velocity Verlet algorithm. Instead, for the Runge-Kutta 4 method, the energy systematically decreases from its exact
value at time 0, meaning that energy conservation is not preserved by the Runge-Kutta 4 method.

We can therefore give the motivation to prefer the velocity Verlet algorithm.

Choice of the velocity Verlet algorithm

The velocity Verlet algorithm is more suitable than the Runge-Kutta 4 method to perform long simulations
of systems with conserved quantities.

Although the velocity Verlet algorithm is symplectic, it can be unstable when trying to solve problems, as shown
in the example below.

We come back to the Cauchy problem studied above:

 2
d x


2
= −ω 2 x,
dt





x(0) = 1,



 dx (0) = 0.



dt
If you solve it using the velocity Verlet algorithm, you get the two coupled equations (with v = dx/dt):

ω 2 h2
 
(h) (h)
xn+1 = xn 1− + hvn(h) ,


 2
ω 2 h2 ω 2 h2
   
 (h) (h) 2 (h)
vn+1 = vn 1− − hω xn 1− ,


2 4
(h)
that you can recast into a simple second-order recursion relation for xn :
(h) (h)
xn+2 = (2 − ω 2 h2 )xn+1 − x(h)
n .

(h)
The solution of this recursion relation is of the form xn = Ar1n + Br2n , with r1 and r2 the two roots of the quadratic
equation
r2 − (2 − ω 2 h2 )r + 1 = 0.
If the discriminant ∆ = (2 − ω 2 h2 )2 − 4 < 0 then the two roots are complex conjugate. This happens if ωh < 2. In
this case, as r1 r2 = 1 (relation between the product of the two roots and the coefficients of the quadratic equation),
they verify |r1 | = |r2 | = 1. You can thus write them as r1 = eiϕ , r2 = e−iϕ . If you finally inject this form into the
quadratic equation, you get
ω 2 h2
cos ϕ = 1 − .
2
(h)
From that, you get that xn = xi cos(nϕ). The method is thus stable.
Instead, if ωh > 2, then the roots r1 and r2 are both real but their product still equals 1. Therefore, one of the roots
(h)
is of absolute value larger than 1 and xn diverges. Therefore, the calculation becomes unstable if h > 2/ω.

We end this section by summarizing the advantages and drawbacks of the velocity Verlet algorithm.

19
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

Velocity Verlet Velocity Verlet

1.5 2.1
Simulation Simulation
1 Exact Exact
2.05
0.5

0 2

E
x

−0.5
1.95
−1

−1.5 1.9
0 2 4 6 8 10 0 20 40 60 80 100
t t

Runge-Kutta 4 Runge-Kutta 4
1.5 2.1
Simulation Simulation
1 Exact Exact
2.05
0.5

0 2
E
x

−0.5
1.95
−1

−1.5 1.9
0 2 4 6 8 10 0 20 40 60 80 100
t t

Figure 4: Energy conservation with the velocity Verlet algorithm and the Runge-Kutta 4 method.
On the left panels, we show the result of the numerical integration of the ODE d2 x/dt2 = −ω 2 x (harmonic
oscillator) obtained with the velocity Verlet algorithm (top) or the Runge-Kutta 4 method (bottom). On the
right panels, we show the energy E = (1/2)m(dx/dt)2 + (1/2)mω 2 x2 as a function of time for both methods on
longer time intervals. While the Runge-Kutta 4 method has a better accuracy than the velocity Verlet method
(the simulation points are close to the exact analytic result), energy systematically decreases although it should
be conserved. Instead, energy oscillates close to the exact value with the velocity Verlet algorithm.

Advantages and Drawbacks of the velocity Verlet method

Advantages:

▶ The velocity Verlet method is relatively simple to implement and only requires to store the values of
position and velocity at the previous step.

▶ The velocity Verlet method is relatively fast (few operations at each time step).

20
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

▶ The velocity Verlet method is more precise than the Euler method but less precise than the Runge-
Kutta 4 method.

▶ The velocity Verlet method is symplectic unlike the Runge-Kutta 4 method: it approximately preserves
conservation laws.

Drawbacks:

▶ The velocity Verlet method can be unstable when you have to solve stiff problems.

VIII. Boundary Value Problems

We end these lecture notes about the numerical integration of ODEs by dealing with the case of Boundary Value
Problems (BVPs), which correspond to ODEs for which all the constants of integration are not given at the same
time. As an illustration, we consider the following problem:
 2 #»
d #»
 
d x #» x


2
= f x, ,
 dt dt



#»
x (ti ) = x#»i , (20)



 #»
x (tf ) = x#»f ,



to be solved in the time interval [ti , tf ]. The problem is that we do not know the initial velocity. Therefore if we
choose it at random, there is no guarantee that at the end of the integration we have #» x (tf ) = x#»f . We propose a
method to determine what initial condition on the velocity we should consider in order to meet the final condition
#»
x (tf ) = x#»f . This method is called the shooting method.

Definition of the shooting method

To solve a BVP given by Eq. (20), we transform it into an optimization problem over IVPs. Said differently,
we apply the following procedure.

1. Integrate the ODE up to time tf with a initial guess for the initial velocity v#»i = d #»
x /dt(ti ).

2. Compare #»
x (tf ) with x#»f .

3. Adjust the initial velocity v#»i and iterate until you meet the condition #»
x (tf ) = x#»f up to a given precision.

For the iteration, you can use a root-finding algorithm applied to the function F : v#»i 7→ #» x (tf ) − x#»f , which
#»
associates to an initial velocity vi the value of the solution at time tf minus the target value. Be careful
that the velocity which allows you to meet the condition #» x (tf ) = x#»f is not necessarily unique (F may have
multiple roots).

In the following, we propose an implementation for a one-dimensional problem with the bisection method as a
root-finding algorithm (see Fig. 5). We assume that we know two values v0 and v0′ of the initial velocity such
that F (v0 ) > 0 and F (v0′ ) < 0, where F : v 7→ x(tf ) − xf with x(t) the numerical solution of the Cauchy
problem
d2 x
 
dx dx
= f x, , x(ti ) = xi , (ti ) = v.
dt2 dt dt

Algorithm for the shooting method with a bisection root-finding scheme

vmin ← v0 ▷ Initiate the left bounding value v0 of the velocity for which F (vmin ) > 0 to v0 .
vmax ← v0′ ▷ Initiate the right bounding value vmax of the velocity for which F (vmax ) < 0 to v0′ .
v ← (vmin + vmax )/2 ▷ Choose the trial velocity as the middle between vmin and vmax .
xend ← final_value_IVP(xi , v) ▷ Compute the value x(tf ) of the solution at time tf for initial con-
ditions x(ti ) = xi and dx/dt(ti ) = v.

21
Benjamin GUISELIN Modélisation et Simulation en Physique (HAP708P)

while |xend /xf − 1| > ϵ do ▷ Iterate as long as the relative difference between x(tf ) and xf is larger than
a given threshold ϵ.
if xend > xf then
vmin ← v ▷ F (v) > 0 so update vmin .
else
vmax ← v ▷ F (v) < 0 so update vmax .
v ← (vmin + vmax )/2 ▷ Choose the trial velocity as the middle between the updated values of vmin
and vmax .
xend ← final_value_IVP(xi , v) ▷ Compute the new value of x(tf ).

F (v) Iteration 1:
v0 + v0′
v1 =
F (v0 ) 2
F (v2 ) Sign change in [v0 , v1 ]
F (v3 )
v1 v0′
v0 v2 v3 v
F (v1 ) Iteration 2:
v0 + v1
F (v0′ ) v2 =
2
Sign change in [v2 , v1 ]

Figure 5: Illustration of the bisection method to find the root of a function. We look for the root of
function F knowing that it lies in the range [v0 , v0′ ] [with F (v0 ) > 0 and F (v0′ ) < 0]. We start with vmin = v0
and vmax = v0′ . At each step we choose v in the middle of the interval [vmin , vmax ] and update vmin to v if
F (v) > 0 or vmax to v if F (v) < 0. After few iterations, v converges to the root of F .

22