Deadline
10:00, Monday 25th November
Task
COMP5930M Scientific Computing
Coursework 2 November 5, 2019
The numbered sections of this document describe problems which are modelled by partial differ ential equations. A numerical model is specified which leads to a nonlinear system of equations. You will use one or more of the algorithms we have covered in the module to produce a numerical solution.
Matlab scripts referred to in this document can be downloaded from Minerva under Learning Resources Coursework. Matlab implementations of some algorithms have been provided as part of the module but you can implement your solutions in any other language if you prefer.
You should submit your answers as a single PDF document via the VLE before the stated deadline. MATLAB code submitted as part of your answers should be included in the document. MATLAB functions should include appropriate help information that describe the purpose and use of that function.
Standard late penalties apply for work submitted after the deadline.
Disclaimer
This is intended as an individual piece of work and, while discussion of the work is encouraged, plagiarism of writing or code in any form is strictly prohibited.
1. A onedimensional PDE: Nonlinear parabolic equation 20 marks total
A nonlinear parabolic equation is written as
u 2u u2
t x2 x 1
for an unknown function ux on spatial domain x 0, 1 and time domain t 0. and are known, positive, constants.
Boundary conditions are specified as u0 0 and u1 1. Initial conditions are specified at t 0 as ux, 0 x.
We numerically approximate 1 using the method of lines on a uniform spatial grid with m nodes on the interval 0, 1 with grid spacing h, and a fixed time step of t.
Code for this problem can be downloaded from Minerva and is in the Q1 folder.
a A numerical approximation using central finite difference methods in space and the
implicit Euler method can be written for the node i as:
uk1uk 2
i i uk1 2 uk1 uk1 uk1 uk1 .
t h2 i1 i i1 4h2 i1 i1
Define the nonlinear system FU 0 that needs to be solved at each time step to obtain a numerical solution of the PDE 1. State the precise number of equations N, and define the unknowns U.
b Solve the numerical problem for N 81 equations on the time interval t 0,2 with 80 time steps. The main program is executed using the runModel.m function. The nonlinear system is implemented in the Matlab function residual.m such that Newtons method can be applied to find a numerical solution. The Jacobian can be computed using fdJacobian.m and linearSolve.m should be used to solve the linear system at each iteration. Note: linearSolve.m uses dense linear algebra.
Plot the solution at t 0 and t 2 and include the plot in your writeup.
How many Newton iterations does the method take in total?
c Solve the same problem for different problem sizes: N 81, 161, 321, 641 choose m
appropriately for each case. Create a table that shows for each different N: i. the total computational time taken, T
ii. the total number of Newton iterations, S
iii. the average time spent per Newton iteration, tS TS
Estimate the algorithmic cost of one Newton iteration as a function of the number of equations N from this simulation data. Does your observation match the theoretical cost of the algorithms you have used?
Hint: Take the values N and the measured tSN and fit a curve of the form
tS CNP by taking logarithms in both sides of the equation to arrive at logtS logCPlogN,
then use polyfit logN, logtS, 1 to find the P that best fits your data.
d Derive the exact expression for row i of the Jacobian for this problem. What is the sparsity pattern of the Jacobian for this problem?
e Repeatthetimingexperimentc,butusingnowtheThomasalgorithmsparseThomas.m to solve the linear system and the tridiagonal implementation of the numerical Ja cobian computation tridiagonalJacobian.m. You will need to modify the code to call newtonAlgorithm.m appropriately.
Create a table that shows for each different N: i. the total computational time taken, T
ii. the total number of Newton iterations, S
iii. the average time spent per Newton iteration, tS TS
Estimate the algorithmic cost of one Newton iteration as a function of the number of equations N from this set of simulations. Does your observation match the theoretical cost of the algorithms you have used?
Hint: Take the values N and the measured tSN and fit a curve of the form
tS CNP by taking logarithms in both sides of the equation to arrive at logtS logCPlogN,
then use polyfit logN, logtS, 1 to find the P that best fits your data.
2. A threedimensional PDE: Nonlinear diffusion 20 marks total
Consider the following PDE for ux, y, z:
2 u 2 u 2 u
x 1u x y 1u y z 1u z gx,y,z, 2 defined on x,y,z 10,103 with ux,y,z 0 on the boundary of the domain, and
some known function gx, y, z that does not depend on u.
Using a finite difference approximation to the PDE 2 on a uniform grid, with grid size h and n nodes in each coordinate direction, the nonlinear equation Fijk 0 to be solved at an internal node may be written as
Fi,j,k
h2
1
1 2
ui,j1,k ui,j,k 2 2
ui,j,k1 ui,j,k 2 2
ui,j1,k ui,j,k ui,j,k1 ui,j,k
1 1
ui,j,k ui,j1,k 2 2
ui,j,k ui,j,k1 2 2
ui,j,k ui1,j,k ui,j,k ui,j1,k
1 ui1,j,k ui,j,k 2
ui1,j,k ui,j,k
ui,j,k ui1,j,k 2 1 2
1
gxi,j,k, yi,j,k, zi,j,k 0
ui,j,k ui,j,k1
Code for this problem can be downloaded from Minerva and is in the Q2 folder. The Matlab function runFDM.m controls the numerical simulation of the discrete nonlinear system 3.
a Identify the finite difference approximations that have been made in the fully discrete form 3. Hint: Use one of the terms, e.g. 1 u2 u , on the lefthandside of
the PDE 2, discretised at a typical internal grid point i,j,k, to illustrate your answer.
3
x x
b Run the model using runFDM.m with grid dimension m 5, 8 and 11 for tol 108 and initial guess ux, 0 0 for all x. Create a table that shows for each different N :
i. the total number of unknowns, N
ii. the total computational time taken, T
iii. the total number of Newton iterations, S
iv. the average time spent per Newton iteration, tS TS
Estimate the algorithmic cost of a single Newton iteration in terms of the size of the nonlinear system N, i.e. T ONp. How does it compare to the theoretical worstcase complexity?
Hint: Take the values N and the measured tSN and fit a curve of the form tS CNP by taking logarithms in both sides of the equation to arrive at
logtS logCPlogN,
then use polyfit logN, logtS, 1 to find the P that best fits your data.
c Replace in runFDM.m the call to linearSolve.m function that uses dense linear alge bra with the alternative matlabLinearSolve.m that uses sparse linear algebra. Solve the problem for m 5 with both dense linear algebra and sparse linear algebra, and show that the numerical solution does not change when we change the solver.
d Replace in runFDM.m the call to a dense Jacobian implementation fnJacobian.m with a call to a sparse implementation, buildJacobian.m. You will need to uncomment the following function calls before calling the Newton algorithm:
A, nz sparsePattern3DFD m ; 3d FDM sparsity pattern
loc, nf sparseJac A,N ; function calls required to evaluate J
i. Visualise the sparsity pattern of A using the command spy in the case m 8 and report the number of nonzero elements nnz
ii. What is the bandwidth of the matrix A in the case m 8? The bandwidth of a matrix A is defined as the smallest number k s.t.
ai,j 0forallijk,
i.e. the smallest distance from the diagonal after which all elements will be 0.
e Use the command lu to compute the LUfactorisation of A in the case m 8. What are the number of nonzero elements L and U, respectively?
Explain what we mean by fillin of the LUfactors.
How could you, in principle, reduce the amount of fillin observed in the factors L and U?
f Solve the problem using the sparse implementation for m 15, 19 and 23. Create a table that shows for each different N:
i. the total number of unknowns, N
ii. the total computational time taken, T
iii. the total number of Newton iterations, S
iv. the average time spent per Newton iteration, tS TS
Estimate the algorithmic cost of this new sparsitybased implementation using the timings of the numerical experiments. Hint: Take the values N and the measured tSN and fit a curve of the form tS CNP by taking logarithms in both sides of the equation to arrive at
logtS logCPlogN,
then use polyfit logN, logtS, 1 to find the P that best fits your data.
Learning objectives
Formulation of sparse nonlinear systems of equations from discretised PDE models. Measuring efficiency of algorithms for large nonlinear systems.
Efficient implementation for sparse nonlinear systems.
Mark scheme
This piece of work is worth 20 of the final module grade.
There are 40 marks in total.
1. 1d PDE 20 marks total
a Formulation of the discrete problem 4 marks
b Simulation 2 marks
c Timings and analysis of complexity, general case 5 marks
d Jacobian structure 4 marks
e Timings and analysis of complexity, tridiagonal case 5 marks
2. 3d PDE 20 marks total
a Formulation of the discrete problem 3 marks
b Simulations with dense version of the solver 5 marks c Comparison of dense and sparse solver outputs 2 marks d Analysis of the sparsity of the Jacobian matrix 2 marks e Analysis of the LU factors of the Jacobian 4 marks
f Simulations with sparse version of the solver 4 marks