5. Higher-Order Partial Differential Equation

For this assignment, you will numerically solve the diffusion equation:

$$\frac{\partial u}{\partial t} =  \frac{\partial^2u}{\partial x^2} $$

Let the \(x\) values be restricted to the interval [0,1], and the boundary values be such that \(u(0,t)=\sin^2(t)\) and \(u(1,t)=0\) (for all \(t\)). The initial value is \(u(x,0)=0\).

The solution should be computed with an explicit time-stepping scheme using timesteps \(\Delta t\) upto a time \(T\) (forward Euler will suffice here), and using a discretization of the interval into \(N\) points \(x_i=i \Delta x\) where \(i=0..N-1\). 

Different from the lecture, you should use the following 4th order approximation to the right hand side of the diffusion equation:

$$\left.\frac{\partial^2u}{\partial x^2}\right|_{x=x_i}\approx\frac{-u_{i-2}+16u_{i-1}-30u_i+16u_{i+1}-u_{i+2}}{12\Delta x^2}$$

where \(i=2..N-3\).

For \(i=1\), this breaks down because \(x_{-1}\) does not exist; use the following instead:

$$\left.\frac{\partial^2u}{\partial x^2}\right|_{x=x_1}\approx\frac{11u_0-20u_1+6u_2+4u_3-u_4}{12\Delta x^2}$$

Similarly, for \(x_{N-2}\), use:

$$\left.\frac{\partial^2u}{\partial x^2}\right|_{x=x_{N-2}}\approx\frac{11u_{N-1}-20u_{N-2}+6u_{N-3}+4u_{N-4}-u_{N-5}}{12\Delta x^2}$$

• For this assignment, you should in fact write three implementations of the simulation:

1. using a full matrix and an explicit loop for the matrix-vector multiplication;
2. using a full matrix and a call to a BLAS routine for the matrix-vector multiplication;
3. using a banded matrix and the appropriate call to a banded matrix-vector BLAS routine.

• You will need to install OpenBLAS in your Teach cluster account.
• During the simulation, the code should make \(P\) snapshots of the value of \(u\) on the \(N\) points to be printed to a file. The format of the output should be such that each line contains the time, the grid spacing, and values of \(u\) at that time (e.g. `std::println(f,"{} {} {}", t, dx, u)` where u is a rvector).
• \(P, N, T, \Delta t\), the name of the output file, and a number 1..3 should be command line parameters of your code. The number should select which of the three implementations will be used. 
• Use the rarray library for your all vectors and matrices. 
• Comment your code.
• Use git
• Use make. 
• The codes should be able to run on the Teach cluster with the gcc/14.3 and rarray/2.8.2 modules loaded.
• Add a 'run' target to your Makefile that runs the case \(P=400\), \(N=200\), \(T=10\), and \(\Delta t=10^{-6}\) (THIS PREVIOUSLY ERRONEOUSLY SAID \(10^{-5}\), WHICH DOES NOT WORK).
• Note that modularity is not a requirement in this assignment, but a few functions to clarify the code structure are expected. 
  

For the banded matrix approach, you have to figure out how to 'pack' the original, 7-banded matrix into a Nx7 full matrix. The following document may be helpful:  https://www.netlib.org/blas/blast-forum/chapter2.pdf.

Submit your repo as a zip file and the output of the 'run' target by March 6, 2026. The usual late policy applies.