Assignment 2 - Solving two 1D problems#

This assignment makes up 20% of the overall marks for the course. The deadline for submitting this assignment is 5pm on Thursday 3 November 2022.

Coursework is to be submitted using the link on Moodle. You should submit a single pdf file containing your code, the output when you run your code, and your answers to any text questions included in the assessment. The easiest ways to create this file are:

  • Write your code and answers in a Jupyter notebook, then select File -> Download as -> PDF via LaTeX (.pdf).

  • Write your code and answers on Google Colab, then select File -> Print, and print it as a pdf.

Tasks you are required to carry out and questions you are required to answer are shown in bold below.

The assignment#

Part 1: Solving a wave problem with sparse matrices#

In this part of the assignment, we want to compute the solution to the following (time-harmonic) wave problem:

\[\begin{split} \begin{align*} \frac{\mathrm{d}^2 u}{\mathrm{d}x^2} + k^2u &= 0&&\text{in }(0, 1),\\ u &= 0&&\text{if }x=0,\\ u &= 1&&\text{if }x=1,\\ \end{align*} \end{split}\]

with wavenumber \(k=29\mathrm{\pi}/2\).

In this part, we will approximately solving this problem using the method of finite differences. We do this by taking an evenly spaced values \(x_0=0, x_1, x_2, ..., x_N=1\) and approximating the value of \(u\) for each value: we will call these approximations \(u_i\). To compute these approximations, we use the approximation

\[ \frac{\mathrm{d}^2u_{i}}{\mathrm{d}x^2} \approx \frac{ u_{i-1}-2u_i+u_{i+1} }{h^2}, \]

where \(h = 1/N\).

With a bit of algebra, we see that the wave problem can be written as

\[ (2-h^2k^2)u_i-u_{i-1}-u_{i+1} = 0 \]

if \(x_i\) is not 0 or 1, and

\[\begin{split} \begin{align*} u_i &= 0 &&\text{if }x_i=0,\\ u_i &= 1 &&\text{if }x_i=1. \end{align*} \end{split}\]

This information can be used to re-write the problem as the matrix-vector problem \(\mathrm{A}\mathbf{u}=\mathbf{f},\) where \(\mathrm{A}\) is a known matrix, \(\mathbf{f}\) is a known vector, and \(\mathbf{u}\) is an unknown vector that we want to compute. The entries of \(\mathbf{f}\) and \(\mathbf{u}\) are given by

\[\begin{split} \begin{align*} \left[\mathbf{u}\right]_i &= u_i,\\ \left[\mathbf{f}\right]_i &= \begin{cases} 1&\text{if }i=N,\\ 0&\text{otherwise}. \end{cases} \end{align*} \end{split}\]

The rows of \(\mathrm{A}\) are given by

\[\begin{split} \left[\mathrm{A}\right]_{i,j} = \begin{cases} 1&\text{if }i=j,\\ 0&\text{otherwise}, \end{cases} \end{split}\]

if \(i=0\) or \(i=N\); and

\[\begin{split} \left[\mathrm{A}\right]_{i, j} = \begin{cases} 2-h^2k^2&\text{if }j=i,\\ -1&\text{if }j=i+1,\\ -1&\text{if }j=i-1.\\ 0&\text{otherwise}, \end{cases} \end{split}\]

otherwise.

Write a Python function that takes \(N\) as an input and returns the matrix \(\mathrm{A}\) and vector \(\mathrm{f}\). You should use an appropriate sparse storage format for the matrix \(\mathrm{A}\).

The function scipy.sparse.linalg.spsolve can be used to solve a sparse matrix-vector problem. Use this to compute the approximate solution for your problem for \(N=10\), \(N=100\), and \(N=1000\). Use matplotlib (or any other plotting library) to plot the solutions for these three values of \(N\).

Briefly (1-2 sentences) comment on your plots: How different are they to each other? Which do you expect to be closest to the actual solution of the wave problem?

This wave problem was carefully chosen so that its exact solution is known: this solution is \(u_\text{exact}(x) = \sin(kx)\). (You can check this by differentiating this twice and substituting, but you do not need to do this part of this assignment.)

A possible approximate measure of the error in your solution can be found by computing

\[ \max_i\left|u_i-u_\text{exact}(x_i)\right|. \]

Compute this error for a range of values for \(N\) of your choice, for the method you wrote above. On axes that both use log scales, plot \(N\) against the error in your solution. You should pick a range of values for \(N\) so that this plot will give you useful information about the methods.

For the same values of \(N\), measure the time taken to compute your approximation for your function. On axes that both use log scales, plot \(N\) against the time taken to compute a solution.

We now want to compute an approximate solution where the measure of error is \(10^{-8}\) or less. By looking at your plots, pick a value of \(N\) that you would expect to give error of \(10^{-8}\) or less. Briefly (1-2 sentences) explain how you picked your value of \(N\) and predict how long the computation will take.

Compute the approximate solution with your value of \(N\). Measure the time taken and the error, and briefly (1-2 sentences) comment on how these compare to your predictions. Your error may turn out to be higher than \(10^{-8}\) for your value of \(N\): if so, you can still get full marks for commenting on why your prediction was not correct. Depending on your implementation and your prediction, a valid conclusion in the section could be “My value of \(N\) is too large for it to be feasible to complete this computation in a reasonable amount of time / without running out of memory”.

Part 2: Solving the heat equation with GPU acceleration#

In this part of the assignment, we want to solve the heat equation

\[\begin{split} \begin{align*} \frac{\mathrm{d}u}{\mathrm{d}t} &= \frac{1}{1000}\frac{\mathrm{d}^2u}{\mathrm{d}x^2}&&\text{for }x\in(0,1),\\ u(x, 0) &= 0,&&\text{if }x\not=0\text{ and }x\not=1\\ u(0,t) &= 10,\\ u(1,t) &= 10. \end{align*} \end{split}\]

This represents a rod that starts at 0 temperature which is heated to a temperature of 10 at both ends.

Again, we will approximately solve this by taking an evenly spaced values \(x_0=0, x_1, x_2, ..., x_N=1\). Additionally, we will take a set of evenly spaced times \(t_0=0,t_1=h, t_2=2h, t_3=3h, ...\), where \(h=1/N\). We will write \(u^{(j)}_{i}\) for the approximate value of \(u\) at point \(x_i\) and time \(t_j\) (ie \(u^{(j)}_{i}\approx u(x_i, t_j)\)).

Approximating both derivatives (similar to what we did in part 1), and doing some algebra, we can rewrite the heat equation as

\[\begin{split} \begin{align*} u^{(j + 1)}_i&=u^{(j)}_i + \frac{u^{(j)}_{i-1}-2u^{(j)}_i+u^{(j)}_{i+1}}{1000h},\\ u^{(0)}_i &= 0,\\ u^{(j)}_{0}&=10,\\ u^{(j)}_{N}&=10. \end{align*} \end{split}\]

This leads us to an iterative method for solving this problem: first, at \(t=0\), we set

\[\begin{split} u^{(0)}_i = \begin{cases} 10 &\text{if }i=0\text{ or }i=N,\\ 0 &\text{otherwise}; \end{cases} \end{split}\]

then for all later values of time, we set

\[\begin{split} u^{(j+1)}_i = \begin{cases} 10 &\text{if }i=0\text{ or }i=N,\\ \displaystyle u^{(j)}_i + \frac{u^{(j)}_{i-1}-2u^{(j)}_i+u^{(j)}_{i+1}}{1000h} &\text{otherwise}. \end{cases} \end{split}\]

Implement this iterative scheme in Python. You should implement this as a function that takes \(N\) as an input.

Using a sensible value of \(N\), plot the temperature of the rod at \(t=1\), \(t=2\) and \(t=10\). Briefly (1-2 sentences) comment on how you picked a value for \(N\).

Use numba.cuda to parallelise your implementation on a GPU. You should think carefully about when data needs to be copied, and be careful not to copy data to/from the GPU when not needed.

Use your code to estimate the time at which the temperature of the midpoint of the rod first exceeds a temperature of 9.8. Briefly (2-3 sentences) describe how you estimated this time. You may choose to use a plot or diagram to aid your description, but it is not essential to include a plot.