Description
AMATH 342 Assignment 1
2. (5 points) Exercise 1.1 of the book by Arieh Iserles. Only prove convergence of the implicit
midpoint rule (1.12). You do not need to prove convergence of the theta method (1.13).
3. Consider the scalar linear problem y
0 = f(y), y(0) = 1, t ∈ [0, t?
], where f(y) = ay with
a < 0.
(a) (1 point) Write down Euler’s method for this problem.
(b) (3 points) Now take a = −2 and t
? = 100. To prove that Euler’s method is convergent,
we proved that
|en,h| ≤ c
λ
(exp(t
?λ) − 1)h,
where λ is the Lipschitz constant associated with f(y) and c a constant bounding the
O(h
2
) term in the Taylor series expansion of y(tn+1) (see the proof of Theorem 1.1 of
the Arieh Iserles book). A reasonable choice is c = 2 (see page 7 of Arieh Iserles book).
Find a value for λ and therefore a bound for the norm of the error of the form |en,h| ≤ αh
with α ∈ R
+. Is this bound of any use in practice?
(c) (4 points) An improved error bound is given by
|en,h| ≤ 1
2
t
?
a
2h.
Prove this error bound. With a = −2 and t
? = 100, compare this error bound to that
found in (b).
Hint: Show that
|en,h| = |(1 + ah)
n − exp(anh)|.
Then use that for −1 x ≤ 0, and n = 0, 1, 2, … that
exp(nx) −
1
2
nx2
exp((n − 1)x) ≤ (1 + x)
n ≤ exp(nx).
2
4. We are given the following ODE:
y
00 + by0 + cy = 0, t ∈ [0, 1], y(0) = 1, y0
(0) = 0,
where b
2 = 4c.
(a) (1 point) Find the exact solution to the above problem.
(b) (1 point) Write the above second-order ODE as a system of first order ODEs.
(c) (5 points) Implement Euler’s method to solve the system of first order ODEs. Take
b = 10. In 3 separate figures plot both the numerical solution and the exact solution. In
Figure 1 use h = 0.5, in Figure 2 use h = 0.05 and in Figure 3 use h = 0.005.
(d) (2 points) Let the error on a grid with time step h be defined as:
Eh = max
0≤nh≤1
|yn,h − y(nh)|
For h = 0.5, h = 0.05 and h = 0.005 compute Eh. You will see that the error satisfies
Eh = O(h
p
). What is p? Explain your answer.
(e) (5 point) Write down the theta-method for the system of first order ODEs of question
(b). Repeat questions (c) and (d) but instead of the Euler method, use the θ-method
with θ = 0.5. Which of the two methods discussed here converges faster to the exact
solution? Explain your answer.
5. The Lorenz equations are a simplified model of convection in the earths atmosphere, and is
given by
dx
dt =σ(y − x)
dy
dt =x(ρ − z) − y
dz
dt =xy − βz
where σ, ρ and β are system parameters.
(a) (5 points) Take σ = 10, ρ = 28 and β = 8/3. Implement Euler’s method to solve
the Lorenz equations. Let t ∈ [0, 50] and take h = 0.002. As initial condition take
(x0, y0, z0) = (1, 0, 0). Plot the solution where the horizontal axis is x and the vertical
axis is z.
(b) (1 point) Repeat question (a) but with ρ = 14.
(c) (2 points) Search online or in a textbook for the definitions of ’chaos’, ’strange attractor’, ’Lorentz attractor’, etc. to explain the plots found in questions (a) and (b). (A
short two or three sentence explanation is sufficient.)
AMATH 342 Assignment 2
2. (a) (5 points) Using Theorem 2.1, derive the three-step Adams-Moulton method.
(b) (5 points) Using Theorem 2.1, derive the three-step Adams-Bashforth method.
3. (a) (5 points) Using Theorem 2.1, derive the two-step Nystrom method.
(b) (5 points) Implement the 2-step Nystrom method to solve the following scalar ODE:
y
0 + y = sin(t
2
), t ≥ 0, y(0) = 0.
Plot y as a function of t for t ∈ [0, 8]. Take h = 0.04.
Regarding the implementation, use Euler’s method to find y1. Use the Nystrom method
for all following time steps.
4. (5 points) Exercise 2.6 of book.
5. (a) (5 points) Construct the Gaussian quadrature formulae for the weight function ω(t) ≡
1, 0 ≤ t ≤ 1, of order two.
(b) (5 points) Construct the Gaussian quadrature formulae for the weight function ω(t) ≡
1, 0 ≤ t ≤ 1, of order four.
(over)
2
6. Consider the linear ODE system
y
0 = Ay where A =
−20 10 0 0
10 −20 10 0
0 10 −20 10
0 0 10 −20
. (1)
Let the initial condition be given by y(0) = [1 1 1 1]T
.
(a) (5 points) Write down the 3-step Adams-Bashforth method for the given linear ODE.
Implement this method to solve (1) over the time interval t ∈ [0, 10]. At each time step
compute the Euclidean norm of the solution and plot this Euclidean norm as a function
of time.
In your implementation, define h = 10/N, with N a positive integer. Given
that the norm of the true solution approaches zero as t → ∞, find out for which values
of N ≥ 50 the numerical method gives you a good approximation.
Regarding the implementation of the 3-step Adams-Bashforth method, use Euler’s method
to find y1 and y2. Use the 3-step Adams-Bashforth method for all following time steps.
(b) (3 points) Repeat part (a) but use a 2-step BDF method instead of the 3-step AdamsBashforth method.
Regarding the implementation of the 2-step BDF method, use the Backward Euler
method to find y1. Use the 2-step BDF method for all following time steps.
(c) (1 point) According to you, which method is best for this problem, the 3-step AdamsBashforth method or the 2-step BDF method? Explain your answer.
7. (5 points) The R¨ossler attractor arose from studying oscillations in chemical reactions. The
equations describing this system are given by
dx
dt = − (y + z)
dy
dt =x + Ay
dz
dt =B + xz − Cz
where A, B and C are system parameters. Take A = 0.2, B = 0.2 and C = 5.7.
Implement a
3-stage ERK method of order 3 to solve the R¨ossler system of equations (you can choose any
3-stage ERK of order 3 that you want, but write down which one you chose). Let t ∈ [0, 500]
and take h = 0.02. As initial condition take (x0, y0, z0) = (0, 0, 0). Plot a 3D figure (with xyz
axes) of the solution.
AMATH 342 Assignment 3
2. (a) (5 points) Exercise 3.8.
(b) (3 points) Explain whether the following three-stage IRK is or is not a collocation method:
k1 = yn + h
1
8
f(tn +
1
4
h, k1) + 1
16 f(tn +
1
2
h, k2) + 1
16 f(tn +
3
4
h, k3)
k2 = yn + h
3
16 f(tn +
1
4
h, k1) + 3
32 f(tn +
1
2
h, k2) + 7
32 f(tn +
3
4
h, k3)
k3 = yn + h
1
2
f(tn +
1
4
h, k1) + 1
8
f(tn +
1
2
h, k2) + 1
8
f(tn +
3
4
h, k3)
yn+1 = yn + h
3
10 f(tn +
1
4
h, k1) + 3
5
f(tn +
1
2
h, k2) + 1
10 f(tn +
3
4
h, k3)
3. (5 points) Exercise 4.4.
4. (5 points) Exercise 4.6a. You do not need to do parts (b) and (c).
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.
AMATH 342 Assignment 4
2. Consider the advection equation:
∂tu + a∂xu = 0, −∞ < x < ∞, t > 0 and u(x, 0) = u0(x), (1)
with a > 0 a constant.
Discretize the x-t plane by choosing a mesh width ∆x and time step
∆t and define the discrete mesh points (xj , tn) by
xj = j∆x, j = …, −1, 0, 1, 2, … and tn = n∆t, n = 0, 1, 2, …
(a) (2 points) Show that the following difference scheme is not stable:
U
n+1
j − U
n
j
∆t
+ a
U
n
j+1 − U
n
j−1
2∆x
= 0. (2)
We will now derive a stable scheme.
(b) (1 point) From (1) it is clear that ∂tu = −a∂xu. Show that ∂ttu = a
2∂xxu.
(c) (2 points) Using the Taylor series expansion of u(x, t + ∆t) around t and question (b),
show that
u(x, t + ∆t) = u(x, t) − ∆ta∂xu(x, t) + 1
2
(∆t)
2
a
2
∂xxu(x, t) + h.o.t. (3)
(h.o.t = higher order terms).
(d) (2 points) Approximating the x-derivatives in (3) by central differences, show that we
obtain the following finite difference scheme:
U
n+1
j = U
n
j −
∆t
2∆x
a
U
n
j+1 − U
n
j−1
+
∆t
2
2∆x
2
a
2
U
n
j+1 − 2U
n
j + U
n
j−1
. (4)
(e) (3 points) Define ν = a∆t/∆x. Using Fourier analysis, find the amplification factor.
Under what condition on ν is the finite difference scheme (4) stable?
(over)
3. Given is the following problem:
∂tu = ∂x (a(x)∂xu), u(0, t) = 0, u(1, t) = 0, u(x, 0) = u0(x), where a(x) > 0. (5)
(a) (2 points) Find the finite difference scheme for (5) using the θ-method in time. You
may use that central differences in space for the ∂x (a(x)∂xu) term is given by
∂x (a(x)∂xu) ≈
aj+1/2
(Uj+1 − Uj ) − aj−1/2
(Uj − Uj−1)
(∆x)
2
. (6)
where aj+1/2 = a(xj+1/2
) and xj+1/2 =
1
2
(xj + xj+1).
(b) (2 points) Show that the maximum principle is satisfied if
2∆t(1 − θ) max
x∈[0,1]
a(x) ≤ (∆x)
2
.
(c) (5 points) Show that (6) is a consistent discretization of ∂x (a(x)∂xu).
4. (4 points) Let 0 < x < 1 and 0 < t < TF and consider the following mixed initial-boundary
value problem:
∂tu = 2∂xxu, u(0, t) = u(1, t) = 1, u(x, 0) = f(x), (7)
where f(x) = 2 if x = 0.5 and f(x) = 1 otherwise.
Suppose the mesh points are chosen to
satisfy
0 = x0 < x1 < x2 < … < xJ−1 < xJ = 1
where xj − xj−1 = ∆x for all j = 1, 2, …, J. Implement the θ-method for problem (7). To
solve the implicit system you need to solve a matrix system of the form AUn+1 = F. Write
down the matrix A and vector F.
Now compute a solution using θ =
1
2
and J = 10.
Make two
plots: in the first plot plot the solution at all computed time levels between 0 and TF = 0.07
using ∆t = ∆x
2
(so plot the solution at t = 0, t = ∆t, t = 2∆t, …, t = TF ); in the second plot,
do the same but use ∆t = 0.5∆x
2
.
Explain the difference in solution.
(over)
3
5. Being able to capture boundary layer effects in fluid dynamics is very important and the grid
plays an important role in being able to capture these effects.
Let 0 < x < 1 and 0 < t < 10 and consider the following advection-diffusion problem:
∂tu + ∂xu −
1
Re
∂xxu = 0, (8)
where Re is the Reynolds number. Consider the following boundary conditions u(0, t) = 1
and u(1, t) = 0 and initial condition u(x, 0) = 1 − x. The exact steady-state solution to this
problem is given by
u(x) = exp(Re) − exp(Re x)
exp(Re) − 1
. (9)
To discretize (8), we first use a uniform grid, where
xi = (i − 1)∆x, i = 1, 2, …, N + 1, (10)
where ∆x = 1/N. As finite difference method, we use, for i = 2, 3, …, N,
U
n+1
i = U
n
i −
∆t
∆x
U
n
i − U
n
i−1
+
1
Re
∆t
∆x
2
U
n
i+1 − 2U
n
i + U
n
i−1
, (11)
(a) (5 points) Take Re = 40. Implement the finite difference discretization (11). For this,
take N = 64 and ∆t small enough such that your scheme is stable. In one figure, plot
the exact steady solution given by (9) and the numerical approximation at time t = 10.
At t = 10, compute the error E = maxi=1,…,N |Ui − u(xi)|. Compute this error also
when N = 128. Based on these two errors, what would you say is the order of the finite
difference method given by (11)?
We will now compute the solution on the following non-uniform grid:
xi =
(
2(1 − c)(i − 1)/N, for i = 1, 2, …, N/2 + 1,
1 − c +
2c(i−1−N/2)
N
, for i = N/2 + 1, N/2 + 2, …, N + 1,
(12)
where c = (2/Re) ln(N).
The reason to use this non-uniform grid is that for high Reynolds
numbers Re the solution rapidly changes in a small neighbourhood of the point x = 1. This
non-uniform grid captures these changes better than a uniform grid.
A finite difference approximation to (8), on the non-uniform grid (12), is given by
U
n+1
i = U
n
i −
∆t
∆xi
U
n
i − U
n
i−1
+
1
Re
∆t
∆xi+1
U
n
i+1 − U
n
i
∆xi+1
−
U
n
i − U
n
i−1
∆xi
, (13)
for i = 2, 3, …, N, and where ∆xi = xi − xi−1.
(b) (3 points) Take Re = 40. Implement the finite difference discretization (13). For this,
take N = 64 and ∆t small enough such that your scheme is stable. In one figure, plot
the exact steady solution given in part (a) and the numerical approximation at time
t = 10.
At t = 10, compute the error E = maxi=1,…,N |Ui − u(xi)|. Compute this error
also when N = 128. Based on these two errors, what would you say is the order of the
finite difference method given by (13)?
(c) (1 point) Now take Re = 1000. In the same figure, plot the solution computed on
the uniform grid and the solution on the non-uniform grid. Use N = 64. Explain the
difference.