1.1 Goal of the project

The goal of this project is to analyze and compare different numerical methods for solving the first order advection PDE.

Following the analytical analysis for stability of the numerical scheme, animation were done to visually illustrate and confirm these results. This was carried for different parameters. The animation was programmed in Mathematica and saved to animated gif files which was then loaded into the HTML version of this report.

Fortran 95 was used for the computation part, while Mathematica was used for the animation and graphics part.

The above link contains all the supporting material for the project, including the Fortran program (in source and windows executable format) used to carry the main computation, and the Mathematica program used to do the animation and the Unix bash file used to process the computation for different parameters.

The specific PDE example used for the analysis and animation was the one provided by Professor Donald Dabdub for the final exam for his MAE 185 course (Numerical methods for mechanical engineers) in spring 2006. This PDE is described below:

Solve numerically

\[ \frac {\partial c}{\partial t} + u \frac {\partial c}{\partial x} = 0 \]

Where \(c(x,t)\) is the concentration of a given material as a function of time and space.

The above is solved for the following 2 cases

1. \(u\) (the advection speed, or the speed at which the mass is being transported) is a constant value given as (2 ft/min.
2. \(u\) is a function of time defined as \(u(t) = \frac {t}{20}\ \text {ft/min}\)

The problem parameters are:

\begin{align*} t & \geq 0\\ 0 & \leq x\leq L \end{align*}

Where \(L=100\ \text {feet.}\). Initial conditions are

\[ c\left ( x,0\right ) =F\left ( x\right ) =\left \{ \begin {array} [c]{ccc}1+\cos \left ( \pi \left ( \frac {x-30}{5}\right ) \right ) & & 25\leq x\leq 25\\ 0 & & \text {otherwise}\end {array} \right . \]

The boundary conditions are

\begin{align*} c\left ( 0,t\right ) & =0\\ c(L,t) & =0 \end{align*}

This PDE is an example of an IBVP (Initial and Boundary Value Problem).

Different numerical methods are used to solve the above PDE. The methods are compared for stability using Von Neumann stability analysis.

The numerical methods are also compared for accuracy. This was done by comparing the numerical solution to the known analytical solution at each time step. The comparison was done by computing the root mean square error (RMSE) between the numerical and the analytical solution at each time step.

The method with the least RMSE at the end of the simulation is considered the most accurate.

The above PDE has a known analytical solution which is

\[ C\left ( x,t\right ) =F\left ( x-ut\right ) \]

The above analytical solution indicates that the initial concentration will move from left to right with the advection speed \(u\).

The formulation of each numerical method is shown below. \(h\) is used to represent \(\Delta x,\) the space between 2 space grid point, or the space step size, and \(\tau \) is used to represent \(\Delta t,\) the time step.

The space line has \(N\) grid points. The spacing \(h\) was fixed at \(0.01\) ft for all the methods and for all the test cases, while \(\tau \) was changed. This made comparing the different methods simpler. The following diagram illustrates the discretization used.

Should we consider the lower left and the lower right grid points above as part of the initial conditions, or part of the boundary conditions?

Stability of each method is derived. Stability is important, since by the Lax-Richtmyer equivalence theorem¹, stability implies convergence of the solution. Convergence of the numerical solutions implies that as the step size becomes smaller, the numerical solution converges to the analytical solution.

Explicit and implicit numerical methods are used. When solving for the future value of the solution at a single node in terms of only past values, the method is called an explicit method. In other words, when the only unknown is the future value of the solution at a single node, and everything else on the right hand side of the finite difference equation is a solution derived at earlier time step, the method is explicit.

An implicit method is one in which the finite difference equation contains the solution at a at future time at more than one node. In other words, future solution are being solved for at more than one node in terms of the solution at earlier time. Implicit methods therefor are usually solved by matrix methods by solving \(Ax=b\) where \(b\) represents present present known solution values, and \(x\) are the unknown future solution values, and \(A\) is the coefficient matrix which will usually be block diagonal (or tri diagonal) in shape.

In the derivations below, the notation of \(C_{i}^{n}\) is used to indicate the solution at time step \(n\) and at space node \(i.\) Hence \(C\left ( x_{i},t_{n}\right ) \) is written as \(C_{i}^{n}\). This notation seems to be more clear than the \(C_{i,n}\) notation.

Different finite difference schemes for solving a PDE are obtained by using different methods of approximating the derivative terms in the PDE.

This will be illustrated using the space derivative \(\frac {\partial c}{\partial x}.\) This derivative can be approximated in one of the following 3 ways (all at time step \(n\))

¹Richtmyer and Morton 1967. p45): "Given a properly posed linear initial value problem and a finite difference approximation to it that satisfies the consistency condition, stability is the necessary and sufficient condition for convergence."