ENGR20005
Numerical Methods in Engineering
Workshop 12
Part A: MATLAB Livescripts
12.1 The livescript ENGR20005 Workshop12p1.mlx runs through the solution of systems of differential equations in MATLAB.
(a) Read through the livescript and make sure you understand what each line of code does.
(b) Modify the livescript to solve the initial value problem dy1 =−expy2
dt
dy2 =−expy1 dt
(12.1a) (12.1b)
with the initial conditions y1(0) = 1 and y2(0) = 2.
12.2 The livescript ENGR20005 Workshop12p2.mlx covers the solution of a common class
of ordinary differential equations called stiff problems.
(a) Read through the livescript and understand how to solve these problems.
(b) Using what you have learnt in this livescript, solve the system of differential equations
(12.2)
=
with the initial conditions y1(0) = 1 and y2(0) = 0.
d y1 dt y2
−100 1 y1 0 1 y2
1
Part B: Applications
12.3 The evolution of N vortices in an inviscid fluid is governed by the following set of differential equations
dxj=−1N ωi(yj−yi) dt 2π r2
(12.3a)
(12.3b)
where (xj,yj) is the position of the jth vortex, |ωj| is the strength of the vortex, the sign of ωj denotes its direction (positive ωj corresponds to an anti-clockwise direction of rotation), and rij is the distance between the ith and jth vortices.
i=1 ij i̸=j
dyj = 1 N ωi(xj − xi) dt2π r2
i=1 ij i̸=j
(xi −xj)2 +(yi −yj)2
(a) Consider a system with 2 counter-clockwise rotating vortices (N = 2) of unit
rij =
strength. Find the differential equations that govern the position of these vortices.
(b) Apply the Explicit Euler method to the system of equations you derived in part (a).
(c) Write a MATLAB function (.m file) that uses your answer to part (b) to solve the system of differential equations you derived in part (a) with the initial conditions (x1(0) = 1,y1(0) = 0) and (x2(0) = −1,y2(0) = 0). Use a time step of ∆t = 0.01 over the interval 0 ≤ t ≤ 200.
(d) Now consider a system with 4 counter-clockwise rotating vortices placed at the points (x, y) ∈ {(1, 1), (1, −1), (−1, −1), (−1, 1)}. Derive the equations of motion from Eq. (12.3).
(e) Apply the RK4 method to the system of differential equations you derived in part (d).
(f) Write a MATLAB function (.m file) that uses the equations you derived in part (e) to solve the system of differential equations. Use a time step of ∆t = 0.01 over the interval 0 ≤ t ≤ 200.
2
12.4 * In the livescript ENGR20005 Workshop9p3.mlx, we saw an example of solving the convection equation
∂φ + U ∂φ = 0 (12.4) ∂t ∂x
with periodic boundary conditions using a backwards difference scheme in space and the explicit Euler method in time.
We’ll now solve Eq. (12.4) with periodic boundary conditions and the initial condition
φ(x, 0) = sin x
over the domain x ∈ [0, 2π]. We’ll also assume an evenly spaced mesh with the spacing ∆x = 2π/(N − 1), where N is the number of nodes.
In order to solve this problem numerically, we need approximate the problem in both space and time.
(a) The first step approximating Eq. (12.4) is to approximate the spatial derivatives. In this case well apply a finite difference approximation to the term ∂
i. Use the backwards difference formula and show that dφi = − U (φi − φi−1)
(12.5)
dt ∆x for the interior nodes i = 1, . . . , N − 1.
ii. Apply the boundary condition to the end nodes and show that the relations at the boundaries are given by
dt ∆x where [C] is the convection matrix given by
dφ0 = − U (φ0 − φN ) dt ∆x
dφN =−U(φN−φN−1) dt ∆x
iii. Hence show that the the discretised problem may be written as d {φ} = U [C]{φ}
(12.6)
−1 0 ··· 0 1 1 −1 .
..
[C]= . ..
0 1−10 0 ··· 0 1 −1
3
∂x
(b) Hopefully you can see that by applying a finite difference approximation to Eq. (12.4), we’ve reduced a partial differential equation to a system of first order ordinary differential equations, which we can then solve using any time stepping method.
i. Apply the explicit Euler method to Eq. (12.6) and show that the iterative formula is given by
{φ}n+1 = {φ}n + U∆t[C]{φ}n (12.7) ∆x
ii. Modify the code snippets given in ENGR20005 Workshop12p1.mlx in order to apply Eq. (12.7). Use your script to solve Eq. (12.4) with a convection velocity of U = 1, a time step of ∆t = 0.01, and 1000 nodes.
12.5 Consider three chemical species A, B, and C which obey the following reaction equa- tions
k1
A 1 −→ A 2
k2
A 2 + A 3 −→ A 1 + A 3
k3
2 A 2 −→ 2 A 3
where k1, k2, and k3 are reaction rate constants.
Using these reaction equations, the time histories of each species’ concentration Ci can
be found by the following set of differential equations. dC1 = −k1C1 + k2C2C3
dt
dC2 = k1C1 − k2C2C3 − k3C2 dt
dC3 =k3C2 dt
We’ll look the case where the reaction constants are
k1 = 0.04, k2 = 10, k3 = 1.5×107
and the initial concentrations are
C1 =0.9, C2 =0.1, C3 =0
(12.8)
(a) Use the explicit Euler method to solve Eq. (12.8) with ∆t = 0.1.
(b) The long-term steady-state species concentration occurs when the change in Ci is
zero. i.e. C ̇i = 0.
Determine the steady-state species concentration by finding the roots of
4
−k1C1 + k2C2C3 0 k1C1 − k2C2C3 − k3C2 = 0
k3C2 0
Hint: Conservation of mass ensures that C1 + C2 + C3 = 1.
Does your answer here match what you found in part (a)? (c) Determine the Jacobi matrix and evaluate it at time t = 0.
(12.9)
(d) Determine the eigenvalues of the Jacobi matrix at time t = 0.
(e) Using your result from part (d), what time step would be required for your answer
in part (a) to be stable.
(f) What time step would be required for stability of the following methods
i. Runge–Kutta order 2 ii. Runge–Kutta order 4
iii. Implicit Euler iv. Crank–Nicolson
v. Heun’s method
vi. Linearised Crank–Nicolson
(g) Verify your findings in part (f) by writing MATLAB functions (.m files) to com- pute the solution.
5