Description
5. To study the trajectories of the stars in the galaxy, H´enon and Heiles developed the following system
of equations
d
2x
dt2
= −
∂V
∂x ,
d
2y
dt2
= −
∂V
∂y , (1)
where V (x, y) is the H´enon-Heiles potential given by
V (x, y) = 1
2
x
2 + y
2
+
x
2
y −
1
3
y
3
.
The energy of the system is defined as E(t) = V (x(t), y(t)) + 1
2
(
dx
dt )
2 + ( dy
dt )
2
. It can be shown that
this energy is conserved during motion.
(a) (1 point) Show that the H´enon-Heiles system (1) can be written as the following first order
system:
dx
dt = p
dp
dt = −x − 2xy
dy
dt = q
dq
dt = −y − (x
2 − y
2
)
(2)
(b) (7 points) Use any 4th order explicit numerical method to find an approximate solution to the
H´enon-Heiles problem (write down which method you use). Consider the following two sets of
initial conditions:
x(0) = 0,
dx
dt (0) = 0.436630946376151, y(0) = 0.095,
dy
dt (0) = 0.03,
and
x(0) = 0,
dx
dt (0) = 0.427001854016272, y(0) = 0.095,
dy
dt (0) = 0.096.
For each set of initial conditions verify that energy is indeed conserved (up to approximately nine
digits behind the comma) on the time interval t ∈ [0, 1000] (as time step, use h = 10−2
). Also,
plot the trajectory in the (x, y) plane for t ∈ [0, 1000].
6. In this question we will solve a model of flame propagation. This problem is described on Cleve Moler’s
blog (Cleve Moler is the inventor of Matlab) and is restated here. If you light a match, the ball of
flame grows rapidly until it reaches a critical size.
Then it remains at that size because the amount
of oxygen being consumed by the combustion in the interior of the ball balances the amount available
through the surface. The model is given by the following non-linear ODE:
y
0 = y
2 − y
3
y(0) = δ t ∈ [0, TF ] , TF =
2
δ
,
where δ describes the initial radius of a ball and y(t) represents the radius of the ball. For very small
δ, this is a stiff ODE. For larger δ, the problem is not very stiff.
(a) (7 points) Implement (explicit) Euler’s method for this problem. Take δ = 0.02. Plot the
solution when h = TF /50 and when h = TF /100. Explain the difference in solution.
(b) (3 points) Take now δ = 0.0001. Plot the solution when h = TF /10000 and h = TF /20000.
Explain what you see.
From (b) we just saw that we need on the order of 20000 time steps for Euler’s method to be stable!
It would be better to use an A-stable method. Below we will implement backward Euler’s method.
Backward Euler’s method for the flame propagation model is given by
yn+1 = yn + hy2
n+1 − hy3
n+1, n = 0, 1, … (3)
We immediately see that we run into a problem here: this is a non-linear equation for yn+1. To solve
this non-linear equation, we use Newton-Raphson.
(c) (3 points) Show that Newton-Raphson applied to (3) results in the following iterative scheme:
w
[i+1] = w
[i] −
F(w
[i]
)
F0(w[i])
, i = 0, 1, …, (4)
where F(w) = w − yn − hw2 + hw3 and F
0
(w) = 1 − 2hw + 3hw2
.
We implement Newton-Raphson as follows:
• First we need a good initial guess. Let’s choose w
[0] = yn.
• We then start iterating over i in (4). We continue to iterate until |w
[i+1] − w
[i]
| < 10−10
.
• Once |w
[i+1] − w
[i]
| < 10−10 is satisfied, we set yn+1 = w
[i+1], and therefore we have solved (3)
up to a tolerance of 10−10, which is usually “good enough”.
Pseudo code for Newton’s method applied to the Backward Euler method for the flame propagation
problem, (3), is given by:
(over)
3
δ = (to be set by user)
TF = 2/δ
h = (to be set by user)
n = 1
t = 0
yn = δ
while t < TF
Set w
[0] = yn
Set ε = 1
while ε > 10−10
Compute F(w
[i]
) = w
[i] − yn − h(w
[i]
)
2 + h(w
[i]
)
3
Compute F
0
(w
[i]
) = 1 − 2hw[i] + 3h(w
[i]
)
2
Compute w
[i+1] = w
[i] − F(w
[i]
)/F0
(w
[i]
)
Compute difference ε = |w
[i+1] − w
[i]
|
Update w
[i] = w
[i+1]
end
Update yn+1 = w
[i+1]
Update t = t + h
Update n = n + 1
end
(d) (10 points) Implement backward Euler’s method as just described. Take δ = 0.0001 and plot the
solution for h = TF /20 and h = TF /100. Given that we expect the transition from y = δ to y = 1
to occur around t = 1/δ, we see that the solutions obtained with h = TF /20 and h = TF /100 are
not very accurate. How small do we need to choose h such that we obtain, at least visually, a
sensible solution?
We learn from this exercise that for stiff ODE’s it may be worth spending a bit more time implementing
an implicit numerical scheme than trying to solve the problem with explicit numerical schemes. We
saw from question (b) that we need on the order of 20000 time steps for Euler’s method to be stable.
For backward Euler’s method we saw from question (c), that even at a time step of h = TF /20 the
method is still stable. Not stability, but accuracy is now the determining factor in choosing our time
step.
4