Franz J. Vesely > CompPhys Tutorial > Partial Differential Equations  
 
 





 
< >
Part II: Ch. 5



5.2 Initial Value Problems II: Conservative-parabolic DE

Diffusion:

$ \frac{\textstyle \partial u}{\textstyle \partial t} = \frac{\textstyle \partial}{\textstyle \partial x} (\lambda \frac{\textstyle \partial u}{\textstyle \partial x}) \;\;\;\; {\rm or} \;\;\; \frac{\textstyle \partial u}{\textstyle \partial t} = \lambda \frac{\textstyle \partial^{2} u}{\textstyle \partial x^{2}} $


Best: 2 second-order schemes

- Crank-Nicholson
- Dufort-Frankel

But: first-order algorithms perform well, too

- FTCS (one up for old Leonhard E.!)
- Implicit (even better - good enough for many purposes!)



Subsections


5.2.1 FTCS Scheme for Parabolic DE

In $ \partial u /\partial t = \lambda \partial^{2} u/\partial x^{2}$, replace $\partial u/\partial t$ by DNGF and $\partial^{2} u/\partial x^{2}$ by DDST formula:
$\Longrightarrow$ "forward time-centered space" algorithm,

$ \begin{eqnarray} \frac{\textstyle 1}{\textstyle \Delta t} \left[ u_{j}^{n+1}-u_{j}^{n}\right] &=& \frac{\textstyle \lambda}{\textstyle (\Delta x)^{2}} \left[u_{j+1}^{n}-2 u_{j}^{n}+u_{j-1}^{n}\right] \end{eqnarray} $

Using $a \equiv \lambda \Delta t / (\Delta x)^{2}$ this is
$ \begin{eqnarray} && u_{j}^{n+1}=(1-2a)u_{j}^{n} + a(u_{j-1}^{n}+u_{j+1}^{n}) \end{eqnarray} $


Stability:

For the $k$-dependent growth factor we find $g(k) = 1-4a\sin^{2}\frac{\textstyle k \Delta x}{\textstyle 2}$, which tells us that for stability the condition is

$ \Delta t \leq \frac{\textstyle (\Delta x)^{2}}{\textstyle 2\lambda} \equiv \tau $

where $\tau$ is the characteristic time for the diffusion over a distance $\Delta x$ (i.e. one lattice space).

The FTCS scheme is simple and stable, but inefficient.


EXERCISE: Remember the thermal conduction problem we considered earlier? If you haven't done it then, do it now, using FTCS. Interpret the behavior of the solution for varying time step sizes in the light of the above stability considerations.




5.2.2 Implicit Scheme of First Order

Take the second spatial derivative at time $t_{n+1}$ instead of $t_{n}$:

$ \frac{\textstyle 1}{\textstyle \Delta t} \left[ u_{j}^{n+1}-u_{j}^{n}\right] = \frac{\textstyle \lambda}{\textstyle (\Delta x)^{2}} \left[ u_{j+1}^{n+1}-2 u_{j}^{n+1}+u_{j-1}^{n+1}\right] $





Again defining $a \equiv \lambda \Delta t / (\Delta x)^{2}$, we find, for each space point $x_{j}$ $\;(j=1,2,..N-1)$,

$ -a u_{j-1}^{n+1}+(1+2a)u_{j}^{n+1}-a u_{j+1}^{n+1}=u_{j}^{n} $


Let the boundary values $u_{0}$ and $u_{N}$ be given; the set of equations may then be written as

$ A \cdot u^{n+1} = u^{n} $

with

$ \begin{eqnarray} A & \equiv & \left( \begin{array}{cccccc} 1 & 0 & 0 & . & . & 0 \\ -a & 1+2a & -a & 0 & . & 0 \\ 0 & . & . & . & 0 & . \\ . & . & . & . & . & . \\ . & . & . & -a & 1+2a & -a \\ . & . & . & 0 & 0 & 1 \end{array} \right) \end{eqnarray} $


$\Longrightarrow$ Solve by Recursion!

Stability: We find

$ g = \frac{\textstyle 1}{\textstyle 1+4a\sin^{2}(k \Delta x/2)} $

Since $\vert g\vert \leq 1$ under all circumstances, we have here an unconditionally stable algorithm!


EXERCISE: Apply the implicit technique to the thermal conduction problem discussed before. Consider the efficiency of the procedure as compared to FTCS. Relate the problem to the Wiener-Levy random walk.




5.2.3 Crank-Nicholson Scheme (CN)

As before, replace $ \partial u/ \partial t $ by $ \Delta_{n} u/\Delta t \equiv (u^{n+1}-u^{n})/\Delta t$.

Noting that this approximation is in fact centered at $t_{n+1/2}$, introduce the same kind of time centering on the right-hand side: taking the mean of $\delta_{j}^{2}u^{n}$ (= FTCS) and $\delta_{j}^{2}u^{n+1}$ (= implicit scheme) we write

$ \frac{\textstyle 1}{\textstyle \Delta t} \left[ u_{j}^{n+1}-u_{j}^{n} \right] = \frac{\textstyle \lambda}{\textstyle 2(\Delta x)^{2}} \left[ (u_{j+1}^{n+1}-2 u_{j}^{n+1}+u_{j-1}^{n+1}) + (u_{j+1}^{n}-2 u_{j}^{n}+u_{j-1}^{n}) \right] $




The Crank-Nicholson formula is of second order in $\Delta t$.

Defining $a \equiv \lambda \Delta t / 2(\Delta x)^{2}$ we may write CN as

$ -a u_{j-1}^{n+1}+(1+2a)u_{j}^{n+1}-a u_{j+1}^{n+1}= a u_{j-1}^{n}+(1-2a)u_{j}^{n}+a u_{j+1}^{n} $


or

$ A \cdot u^{n+1} = B \cdot u^{n} $

with

$ A \! \equiv \! \left( \begin{array}{cccccc} 1 & 0 & 0 & . & . & 0 \\ -a & 1\!+\!2a & -a & 0 & . & 0 \\ 0 & . & . & . & 0 & . \\ . & . & . & . & . & . \\ . & . & . & -a & 1\!+\!2a & -a \\ . & . & . & 0 & 0 & 1 \end{array} \right) \;\; B \! \equiv \! \left( \begin{array}{cccccc} 1 & 0 & 0 & . & . & 0 \\ a & 1\!-\!2a & a & 0 & . & 0 \\ 0 & . & . & . & 0 & . \\ . & . & . & . & . & . \\ . & . & . & a & 1\!-\!2a & a \\ . & . & . & 0 & 0 & 1 \end{array} \right) $



Tridiagonal $\Longrightarrow$ Solve by Recursion!

Stability of CN:

The amplification factor is
$ g(k)=\frac{\textstyle 1-2a\sin^{2}(k\Delta x/2)}{\textstyle 1+2a\sin^{2}(k\Delta x/2)} \leq 1 \; , $

which makes the CN method unconditionally stable.


EXAMPLE: The time-dependent Schroedinger equation,

$ \frac{\textstyle \partial u}{\textstyle \partial t} = -i H u \; , \;\;\; {\rm with} \;\; H \equiv \frac{\textstyle \partial^{2}}{\textstyle \partial x^{2}} + U(x) $

when rewritten à la Crank-Nicholson, reads

$ \frac{\textstyle 1}{\textstyle \Delta t}[u_{j}^{n+1}-u_{j}^{n}] = - \frac{\textstyle i}{\textstyle 2} [(Hu)_{j}^{n+1}+(Hu)_{j}^{n}] = - \frac{\textstyle i}{\textstyle 2} \left[ \frac{\textstyle \delta_{j}^{2}u_{j}^{n+1}}{\textstyle (\Delta x)^{2}} +U_{j}u_{j}^{n+1} +\frac{\textstyle \delta_{j}^{2}u_{j}^{n}}{\textstyle (\Delta x)^{2}} +U_{j}u_{j}^{n} \right] $

With $a \equiv \Delta t/2 (\Delta x)^{2}$ and $b_{j} \equiv U(x_{j}) \Delta t/2$ this leads to

$ \begin{eqnarray} (ia) u_{j-1}^{n+1}+(1-2ia+ib_{j})u_{j}^{n+1}&\!+&\!(ia) u_{j+1}^{n+1}= \\ && =(-ia) u_{j-1}^{n}+(1+2ia-ib_{j})u_{j}^{n}+(-ia) u_{j+1}^{n} \end{eqnarray} $

Again, we have a tridiagonal system which may be inverted very efficiently by recursion.



5.2.4 Dufort-Frankel Scheme (DF)

DST in time and space, but in place of $-2u_{j}^{n}$ use $-(u_{j}^{n+1}+u_{j}^{n-1})$:

$ \frac{\textstyle 1}{\textstyle 2 \Delta t} \left[ u_{j}^{n+1}-u_{j}^{n-1} \right] = \frac{\textstyle \lambda}{\textstyle (\Delta x)^{2}} \left[ u_{j+1}^{n}-(u_{j}^{n+1}+u_{j}^{n-1})+u_{j-1}^{n} \right] $

or, with $a \equiv 2 \lambda \Delta t/(\Delta x)^{2}$,

$ \begin{eqnarray} && u_{j}^{n+1} = \frac{\textstyle 1-a}{\textstyle 1+a}u_{j}^{n-1} + \frac{\textstyle a}{\textstyle 1+a} \left[ u_{j+1}^{n}+ u_{j-1}^{n}\right] \end{eqnarray} $




The DF algorithm is of second order in $\Delta t$. In contrast to CN, it is an explicit expression for $u_{j}^{n+1}$.


Stability:

$ g =\frac{\textstyle 1}{\textstyle 1+a} \left[ a \cos k \Delta x \pm \sqrt{1-a^{2}\sin^{2}k \Delta x} \right] $

Considering in turn the cases $a^{2}\sin^{2}k \Delta x \geq 1$ and $\dots < 1$ we find that $|g|^{2} \leq 1$ always; the method is unconditionally stable.



5.2.5 Resumé: Conservative-parabolic DE


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


- Use Crank-Nicholson! (2nd order implicit scheme; needs recursion)

- If too lazy for implicit scheme, use Dufort-Frankel (2nd order explicit)

- If 1st order is sufficient, use implicit scheme


vesely 2006-01-23

< >