ENGR20005
Numerical Methods in Engineering
Workshop 9
Part A: MATLAB Livescripts
9.1 The livescript ENGR20005 Workshop9p1.mlx runs through the solution of boundary value problems in MATLAB.
(a) Read through the livescript and make sure you understand what each line of code does.
(9.1)
(b) Modify the livescript to solve the the boundary value problem d2y = 4(y − x)
dx2
with the boundary conditions y(0) = 0 and y(1) = 2.
Compare your solution with the analytical solution
y(x) = e2(e4 − 1)−1(e2x − e−2x) + x
9.2 The livescript ENGR20005 Workshop9p2.mlx runs through the solution of boundary value problems using finite difference methods.
(a) Read through the livescript and make sure you understand what each line of code does.
(b) Modify the livescript to solve Eq. (9.1).
9.3 The livescript ENGR20005 Workshop9p3.mlx runs through the different types of bound- ary conditions commonly used in the solution of boundary value problems. Read through the livescript and make sure you understand how to implement each type.
9.4 The livescript ENGR20005 Workshop9p4.mlx runs through the solution of boundary value problems using spectral collocation methods.
(a) Read through the livescript and make sure you understand what each line of code does.
(b) Modify the livescript to solve Eq. (9.1). Compare your answer with the finite difference method in Question 9.2.
1
Part B: Problems
9.5 Consider the following differential equation
d2y + y = 0
dx2
(a) By pairing Eq. (9.2) with the following boundary conditions i. y(0) = 1 and y(π/2) = 1.
ii. y(0)=1andy(π)=1. iii. y(0) = 1 and y(2π) = 1.
(9.2)
solve the boundary value problem using the finite difference method. Use as many points as you think is necessary.
(b) Repeat part (a) using the spectral collocation method.
(c) Solve Eq. (9.2) analytically with the boundary conditions listed in part (a). Compare the analytical solution with your answers in part (a).
9.6 The flow of fluid between two plates, where one of them is stationary and the other is moving with velocity U, is given by the boundary value problem
d2u + α = 0 (9.3) dy2
with the boundary conditions u(0) = 0 and u(1) = 1.
NOTE: The dimensionless parameter α = Gd2 measures the relative strength between
pressure and viscous forces.
(a) Solve Eq. (9.3) analytically.
μU
(b) Assuming that α = 1, discretise Eq. (9.3) using the central difference scheme with 5 evenly spaced nodes.
(c) Solve the set equations you obtained in part (b). Compare your solution with your answer in part (a).
(d) Repeat the problem with the spectral collocation method with 5 Gauss–Lobatto nodes. Compare your answer with the finite difference method.
(e) Repeat the problem with α = 0.01 and 100.
2
9.7 * In this question, you will learn the steps to be taken to calculate partial double derivatives of a two-dimensional function. You will learn how you can extend the derivative matrices idea from one-dimension to calculate partial derivatives of two- dimensional functions. We will first construct the matrix for one-dimensional function and then extend the analysis to two-dimensions. In this example, we will calculate the double partial derivatives. The calculation of partial derivatives follows very similar steps.
As was shown in lectures, the double derivative of a function uxx(x) of a function u(x) can be calculated using the second order central, forward and backward difference schemes
uxx(xi) = (ui−1 − 2ui + ui+1)/∆2x uxx(x0) = (u0 − 2u1 + u2)/∆2x uxx(xNx)=(uNx −2uNx−1+uNx−2)/∆2x.
u(x) is evaluated at equally spaced grid points x = x0, x1, ..xNx, i.e. xi − xi−1 = ∆x. ui = u(xi) and i = 0..Nx where Nx + 1 is the number of grid points in x. The above equations can be written more concisely in matrix form as {uxx} = [Dx2]{u}. The explicit structure for the matrix equations for Nx = 4 and Nx = 8 are shown below
{uxx} =
uxx(x4) = 0 0 ∆2
uxx(x5) x 0 0
0 1−21 0 0 0 1 −2 1
(9.5)
uxx(x0) 1 −2 1 0 0 u0
1 uxx(x1) 1 −2 1 0 0 u1
{uxx}= uxx(x2) = 0 1 −2 1 0 u2 ∆2 uxx(x3) x 0 0 1 −2 1u3
(9.4)
uxx (x4 ) 0 0 1
− 2
1
u 4
uxx(x0) 1 −2 1 0 0
uxx(x1) 1 −2 1 0 0
0 1 −2 1 0 0 0 1−21
0 0 0 0 0 0
0 0 u1
0 0 u2
0 0 u3
0 0 u4
uxx(x2)
1 uxx (x3 )
0 0 u5
xx6 6 u (x) 0 0 0 0 0 1 −2 1 0 u
xx7
u (x) 0 0 0
uxx (x8 ) 0 0 0
0 0 0 1 −2 0 0 0 1 −2
7 1 u
1 u 8
(a) Write a MATLAB program to compute the [Dx2] matrix and use it to calculate and plot the double derivative of u(x) = x7 by performing the matrix multiplication, {uxx}=[Dx2]×{u},inMATLAB.Usethedomainx∈[−1,1]. Checkthatyour answer is correct by comparing with the exact solution uxx(x) = 42×5. How many grid points do you need to use to get a good approximation for uxx(x)? Note that the one sided approximation of uxx(x) is not very good at the end points.
3
0
0
0 0 u0
(b) Toseehow[Dx2]canbeusedtocalculatethedoublederivativeofatwo-dimensional function uji = u(xi,yj), we will write the two-dimensional matrix uj,i in vector form
{u} = (u0,u1,…,uNx,uNx+1,ul,…,uNx×(Ny+1))T = (u0,0,u0,1,…,u0,Nx,u1,0,ui,j,…,uNy,Nx)T
where l = i + jNx, i = 0..Nx and j = 0..Ny. Look at section L16.1 of the lecture slides if you are still unsure of how the data is organised.
(c) Assume that the mesh has 5×5 grid points. Show that the two-dimensional double partial derivative matrices in the x and y directions, [Dx2] and Dy2 are given by Eqs. (9.15) and (9.16) (see the last two pages of this document) respectively.
The result above can be extended to calculate the partial double derivatives on a Nx × Ny mesh
uxx(xi, yi) = [I] ⊗ [D2]{u} (9.6) uyy(xi, yi) = [D2] ⊗ [I]{u} (9.7)
where ⊗ is some operation on the one-dimensional second derivative matrix [D2] and the identity matrix [I] that produces the correct two-dimensional derivative matrix. The operation ⊗ on two matrices is called the Kronecker product which is implemented in MATLAB using the kron() function. For an m × n matrix [A] and a p×q matrix [B] is defined as
[A] ⊗ [B] = . . an1[B] · · ·
(d) Use your matrices above to calculate
uxx (xj , yi )
and
. ann[B].
(9.8)
(9.9)
a11[B] · · · a1n[B] . .. .
uyy(xj,yi) (9.10) where u(x, y) = sin(y)x5 for x ∈ [−1, 1] and y ∈ [0, 2]. Plot your computed uxx
and uyy and compare with the exact solution uxx = 20 sin(y)x3
4
and
uyy =−sin(y)x5
(e) Combine the two operations above to calculate and plot
∂2u + ∂2u ∂x2 ∂y2
(9.11)
where u(x,y) = sin(y)x5 for x ∈ [−1,1] and y ∈ [0,2].
9.8 An idealised model of the displacement of a vibrating drum is given by the wave
equation
∂2u = c2∂2u + ∂2u (9.12) ∂t2 ∂x2 ∂y2
Using some mathematical magic, part of the solution of Eq. (9.12) may be found by solving the Helmholtz equation
−∂2φ + ∂2φ = k2φ (9.13) ∂x2 ∂y2
where k is a constant.
Consider an L-shaped domain with the following grid points
(a) Use the central difference method to discretise Eq. (9.13) on the L-shaped domain above. Assume that φ = 0 along the boundary.
(b) Write the system of equations you found in part (a) in the form
−[L]{φ} = k2[I]{φ} (9.14)
5
Notice how this is an eigenvalue problem.
(c) Use eig() to determine the eigenvalues and eigenvectors of [L].
These are the vibrational modes of an L-shaped drum, which are essentially the sounds that such a drum can make.
(d) Plot the first 3 eigenvectors of Eq. (9.14).
(e) It turns out the MATLAB logo is given by the eigenvector associated with the first eigenvalue. Repeat the problem with additional nodes and reproduce the MATLAB logo.
6
7
(9.15)
1 −2 1 0 0 0 0 0 1 −2 1 0 0 0 0 0 0 1 −2 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0000 0 0000 0 0000 0 0000 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0000 0 0000 0 0000 0 0000 0 0 0 0 0 0 0000 0 0000 0 0000 0 0000 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0000 0 0000 0 0000 0 0000 0 0 0 0 0 0 0000 0 0000 0 0000 0 0000 0 0000 0 0 0 0 0 0 0000 0 0000 0 0000 0
0 0 0 0 0 0 00 00 00 00 0 0 00 00 00 00
00
00
00 0
00 0
0 0 0
00 0
00 0
00 0
00 0
001−2 1
001−2 1
0 0 0 1 −2 1 0 0 0 0
[D2]=1 00 0 x∆2x
0000 0000 0000 0 0000 0 0000 0 0 0 0 0 0 0000 0 0000 0 0000 0 0 0 0 0 0 0000 0 0000 0 0000 0 0000 0 0 0 0 0 0 0000 0 0 0 0 0 0
0001−21000 0
00
0 0 0
00 0
00 0
00 0
0 0 0
00 0
00 0
00 0 00 0 0 0 0 00 0
0 0 0 0 0000 0000 0000 0 0 0 0 0000 0000 0000 0000 0 0 0 0 0000 0 0 0 0
1 −2 1 0 0 0
1−2100 0
0 001−2 1
0 001−2 1
0 0 0 0 1 −2 1 0 0 0 0
0 0 00 00 00 0 0 00 00 00 00
0 0 0
0 0000
0 0000
0 0000 0 0 0000 0 0 0 0 0 0 0 0 0000 0 0 0 0 0 0 0
0 0 0 1 −2 1 0 00001−21 0 0 0 0 1 −2 1
1−2100 0 1−2100 0
1−2100 0 1−2100 0 001−2 1 001−2 1
1−2100 0 1−2100 0 001−2 1 001−2 1
8
(9.16)
[D2] = 1 y∆2y
0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0
1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 − 2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1
0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 −2 0 0 0 0 1 0 0 0 0 0 0 0 0