ENGR20005
Numerical Methods in Engineering
Workshop 11
Part A: MATLAB Livescripts
11.1 The livescript ENGR20005 Workshop11p1.mlx runs through the solution of initial value problems using implicit methods.
(a) (b)
Read through the livescript and make sure you understand what each line of code does.
Modify the livescript to solve the the initial value problem
dx = − exp x (11.1)
dt with the initial conditions x(0) = 1.
11.2 The
value problems using Runge–Kutta methods.
livescript ENGR20005 Workshop11p2.mlx runs through the solution of initial (a) Read through the livescript and make sure you understand what each line of code
does.
(b) Modify the livescript to solve Eq. (11.1)
1
Part B: Applications
11.3 In Workshop 10, we looked at the logistic equation for modelling population growth in a finite community
(a)
dN =rN1−N (11.2) dt K
Using the code snippets in ENGR20005 Workshop11p2.mlx, apply the RK2 method to solve Eq. (11.2). Assume that the population is initially N(0) = 10, the growth rate is r = 2, and the carrying capacity is K = 1000.
(b) i.
To apply the Crank–Nicolson method to Eq. (11.3), show that at time step tn we need to solve the nonlinear equation
N =N +∆trN 1−Nn+rN 1−Nn+1 (11.3) n+1 n 2 n K n+1 K
ii. Apply the Newton–Raphson method to Eq. (11.3) and show that at each step Nn+1 = p is given by the iterative formula
(11.4)
(11.5)
p(l) − ∆trp(l)1− p(l) −N − ∆trN 1− Nn (l+1) (l) 2 K n 2 n K
p = p − ∆t 2rp(l) 1−2r−K
iii. Modify you script from part (a) to apply Eq. (11.4). 11.4 * Consider the initial value problem
with the initial condition x(0) = 0.
(a)
(b) (c)
Show that the analytical solution is given by
x(t) = 4t5/4 5
dx = x1/5 dt
Use the explicit Euler method with a time step of ∆t = 0.01 to solve Eq. (11.5) and show that it fails to find the correct solution. Why is this the case?
i. Apply the implicit Euler method to Eq. (11.5) and show that at each step, you have the root finding problem
where p = xn+1.
h(p)=p−∆tp1/5 −xn =0 (11.6)
2
ii. Apply the Newton–Raphson method to Eq. (11.6). Hence argue that if we choose the initial guess p(0) = xn, the method will fail to converge to the correct root.
iii. Instead of using Newton–Raphson to find the root of Eq. (11.6), we’ll now consider using fixed point iteration.
A. Show that possible fixed point formulae for Eq. (11.6) are I. g(p) = xn + ∆tp1/5
II. g(p) = xn 1−∆tp−4/5
B. To see which fixed point formula we should use, consider the initial time step n = 0 so that x0 = 0. We’ll also consider the same step size as in part (b) ∆t = 0.01. Apply your knowledge from the livescript ENGR20005 Workshop2p2.mlx to show that only
g(p) = xn + ∆tp1/5
will converge to the correct value.
Hint: The solutions to h(p) = 0 in this case are p = 0 and p = 1/(100 10). iv. Modify the code snippets on ENGR20005 Workshop11p1.mlx to solve Eq.
(11.5).
(d) Using your answer to part (c), how can you modify your script in part (b) to provide an accurate solution to Eq. (11.5).
11.5 The Explicit Euler, Implicit Euler, and Crank–Nicolson methods are special cases of the very creatively named theta method
xn+1 =xn +∆t[θf(xn,tn)+(1−θ)f(xn+1,tn+1)] (11.7)
where 0 ≤ θ ≤ 1.
Note that we recover the explicit Euler method for θ = 1, the Crank–Nicolson method
for θ = 1/2, and the implicit Euler method for θ = 0.
(a) Determine the truncation error of the theta method and use this result to deter-
mine the order of the method for 0 ≤ θ ≤ 1.
(b) Find an expression for the region of absolute stability for arbitrary θ.
11.6 One of the main difficulties with the Crank–Nicolson method
xn+1 =xn+∆t[f(xn,tn)+f(xn+1,tn+1)] (11.8) 2
is that it is an implicit method, meaning that for nonlinear ODEs, we have to solve a set of nonlinear equations.
3
√
One possible way of making Eq. (11.8) explicit is to use the explicit Euler method to guess the value of f at tn+1, giving Heun’s method
(11.9)
(b) Determine the region of absolute stability of Eq. (11.9).
(c) Compare Eq. (11.9) with the general RK2 method and determine the constants
a1, a2, p1, and q11.
A second way of making the Crank-Nicolson method explicit is to linearise it, giving
(11.10)
x ̃n+1 =xn +∆tf(tn,xn)
xn+1 =xn+∆t[f(tn,xn)+f(tn+1,x ̃n+1)]
2
(a) Determine the order of accuracy of Eq. (11.9).
xn+1 = xn + 2
(d) Determine the order of accuracy of Eq. (11.10).
(e) Determine the region of absolute stability of Eq. (11.10).
∆t(f(xn,tn+1)+f(xn,tn))
4
1 − ∆t ∂f |(xn,tn+1) 2 ∂x