PROJECT I

Write a computer code to solve the one-dimensional scalar diffusion equation:

$$\frac{\partial }{\partial t} A(x,t) = \kappa \frac{\partial ^2}{\partial x^2} A(x,t),$$ (5.1)

with zero boundary conditions at $x = \pm\infty$.

Use the explicit FTCS method with Dirichlet boundary conditions, of form $A(-L/2,t) = A(L/2,t) = 0$, and the initial condition of $A(x,0) = \frac{4}{\sqrt{2\pi}}\exp(-8 x^2)$. Note that the `window' of $x = \pm L$ must be chosen to be large enough so that having a finite window doesn't introduce a large error, but not too large so that the step-size in x produces large discretisation errors.

The code should be written in Matlab; as an option for experienced students, you may wish to convert this to C, to demonstrate the speed improvement obtained with compiled code. An example code in Matlab is attached to these notes. It will need some modifications to the initial conditions! You should try to rewrite this code in a matrix equation form.

The C code you write can either plot the data with PLPLOT, or you can read the output into a Matlab program for plotting. You could also even try using the C-subroutine options of Matlab. The purpose of this is to show how Matlab can be used for program development and prototyping, with C code generally giving faster production runs after the prototyping stage.

The sample code deliberately plots the results for only a subset of all the calculated time-points--do you understand why this is? Note that the analytic solution to this problem (on an infinite line) is something like (check this yourself):

$$A(x,t) = \frac{1}{\sigma (t) \sqrt{(2\pi)}} \exp\biggl [\frac{-x^2}{2 \sigma^2(t)} \biggr],$$ (5.2)

where $\sigma (t) = \sqrt{2\kappa t + 1/16}$.

Your report should include:

• A general description of the problem, and the physics background.
• An FTCS code that you wrote and tested, with a graph of results and errors.
• A table showing how maximum errors vary with step-size and window-size, including what happens when $\Delta x^2/\Delta t$ is varied.
• Results for an infinite domain, using the implicit FTCS method. This is very simple if you use the Matlab matrix notation,
• The source code should also be attached to an email to cochrane@physics.uq.edu.au, along with the relevant input parameters used to generate the plots.