COMP5930M - Scientific Computing Coursework 2
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
COMP5930M - Scientific Computing
Coursework 2
September 21, 2022
Deadline
09:00, Monday 19th December
Task
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 Minerva 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 one-dimensional PDE: Nonlinear parabolic equation
[10 marks total]
Consider the nonlinear parabolic PDE: find u(x, t) such that
= i + α ╱ 、2
(1)
in the spatial interval x ∈ (0, 1) and time domain t > 0.
Here i and α are known, positive, constants.
Boundary conditions are specified as u(0) = 0 and u(1) = 1.
Initial conditions are specified at t = 0 as u(x, 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) Find the fully discrete formulation for (1) using the central finite difference formulas
∂u ui+1 - ui − 1 ∂2u ui+1 - 2ui + ui − 1
Define the nonlinear system F(U) = 0 that needs to be solved at each time step to obtain a numerical solution of the PDE (1).
(b) Derive the exact expression for each of the non-zero elements in row i of the Jacobian for this problem.
(c) Solve the problem for different problem sizes: N = 81, 161, 321, 641, where N is the number of unknowns to solve for (choose the grid size m appropriately to get the correct N). Use NT = 40 time steps to solve from t0 = 0 to tf = 2. Create a table that shows the total computational time taken T, the total number of Newton iterations S, and the average number of Newton iterations per time step TS = T/S. Time the execution of your code using Matlab functions tic and toc. Show 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 = T/S
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 tS (N) and fit a curve of the form
tS = CNP by taking logarithms in both sides of the equation to arrive at log tS = log C + P log N,
then use polyfit( log(N), log(tS), 1 ) to find the P that best fits your data.
(d) Repeat the timing experiment (c), but using now the Thomas algorithm sparseThomas .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 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 = T/S
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 tS (N) and fit a curve of the form
tS = CNP by taking logarithms in both sides of the equation to arrive at
log tS = log C + P log N,
then use polyfit( log(N), log(tS), 1 ) to find the P that best fits your data.
2. A three-dimensional PDE: Nonlinear diffusion
[10 marks total]
Consider the following PDE for u(x, y, z):
╱ (1 + u2 ) 、 + ╱ (1 + u2 ) 、 + ╱ (1 + u2 ) 、 = g(x, y, z), (2)
defined on [x, y, z] ∈ [-10, 10]3 with u(x, y, z) = 0 on the boundary of the domain, and some known function g(x, 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 = ┌╱ 1 + ╱ 、2 \ (ui+1,j,k - ui,j,k) - ╱ 1 + ╱ 、2 \ (ui,j,k - ui − 1,j,k) + ╱ 1 + ╱ 、2 \ (ui,j+1,k - ui,j,k) - ╱ 1 + ╱ 、2 \ (ui,j,k - ui,j − 1,k) + ╱ 1 + ╱ 、2 \ (ui,j,k+1 - ui,j,k) - ╱ 1 + ╱ 、2 \ (ui,j,k - ui,j,k − 1 )┐
- g(xi,j,k, yi,j,k, zi,j,k) = 0 (3)
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). To solve the problem call runFDM(m), where m is the number of grid points in each direction so that the total number of unknowns is N = m3 .
(a) The Jacobian matrix is computed in the routine newtonAlgorithm. Extract the Jacobian matrix J(x0 ) from the first Newton iteration of the algorithm. Then:
i. Visualise the sparsity pattern of the Jacobian matrix using the command spy in the case m = 10 and include the plot in your report. How many non-zero elements does it have?
ii. Use the command lu to compute the LU-factorisation of A in the case m = 10. Visualise the sparsity pattern of the LU-factors using the command spy and include the plot in your report. What are the number of non-zero elements L and U , respectively?
iii. Use these sparsity patterns to explain what we mean by hll-in of the LU-factors. Do you observe fill-in here?
(b) Solve the tangent problem by using the iterative linear solver iterativeLinearSolve .m based on the GMRES iteration (implemented by the gmres command in MATLAB).
Use a fixed linear tolerance of linTol = 10 −7 and 100 maximum linear iterations. Solve the problem for m = 10, 15, 20 and measure in each case the number of Newton iterations and total number of linear iterations until convergence is achieved. Perform the same experiment with two different GMRES strategies:
i. No preconditioner (default);
ii. Inexact LU-preconditioner M ≈ LU computed by the MATLAB command
options .type = ’crout’;
options .milu = ’row’;
options .droptol = 0 .05;
[L,U] = ilu(A,options);
M=L*U;
(modify the file iterativeLinearSolve .m to implement the preconditioner).
Create a table that measures the total number of Newton iterations and total number of linear iterations in each case for m = 10, 15, 20 for both choices of preconditioning strategies. How does preconditioning affect the number of linear iterations required?
(c) Modify the file iterativeLinearSolve .m to implement an inexact Newton-Krylov
method where the GMRES linear stopping tolerance is chosen as
linTol = maxeη )F(xk)), 10 −7 ( ,
where F(xk) is the nonlinear residual at iteration k and the parameter η is chosen from the Eisenstat-Walker rule:
η = 0.1 )F(xk))2
Repeat the experiment from (b) and update the table to include the number of Newton iterations and total number of linear iterations in this case. How does the inexact Newton-Krylov influence the total number of iterations required in this case?
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 20 marks in total.
1. One-dimensional PDE [10 marks total]
(a) Formulation of the discrete problem [3 marks]
(b) Jacobian structure [3 marks]
(c) Timings and analysis of complexity, general case [2 marks]
(d) Timings and analysis of complexity, tridiagonal case [2 marks]
2. Three-dimensional PDE [10 marks total]
(a) Analysis of the sparsity of the Jacobian matrix and the the LU factors [4 marks]
(b) Experiments with the iterative GMRES linear solver [2 marks]
(c) Experiments with the inexact Newton-Krylov method [4 marks]
2022-12-15