Higher Order Finite Difference Schemes (Part 4)

by

Privious part: Higher Order Finite Difference Schemes for the Heat Eq (Part 3)

Higher Order Finite Difference Schemes for the Heat Eq (Part 2)

Higher Order Finite Difference Schemes for the Heat Eq (Part 1)

Schemes using Three Points in Time and Three Points in Space

Higher Order Finite Difference Schemes for the Heat Eq (Part 4)

The numerical solution for the equation (??) is calculated by using FTCS, BTCS, Crank Nicolson and Adams Bashforth schemes using three points formulae, which are given in details.

FTCS

The approximation of time derivative in (??) complimented with forward difference

(1)   \begin{equation*} \dfrac{\partial w}{\partial t}\big|{t{m},x_i}=\dfrac{-w_i^{m+2}+4w_i^{m+1}-3w_i^m}{2\Delta t}+\mathcal{O}(\Delta t^2)\end{equation*}


The approximation Central difference is used to approximate
(\dfrac{\partial^2 w}{\partial x^2}){x_i} and calculate all the terms for a supposed time m

(2)   \begin{equation*} \dfrac{\partial^2 w}{\partial x^2}\big|{x_i} =\dfrac{w_{i+1}^m-2w_i^m+w_{i-1}^m}{\Delta x^2}+\mathcal{O}(\Delta x^2)\end{equation*}


Using these equations (1) and (2) in (??) and combine the terms of truncation error to obtain

(3)   \begin{equation*} \dfrac{-w_i^{m+2}+4w_i^{m+1}-3w_i^m}{2\Delta t}=\alpha \big[\dfrac{w_{i+1}^m-2w_i^m+w_{i-1}^m}{\Delta x^2} \big]+ \mathcal{O}(\Delta t^2 , \Delta x^2)\end{equation*}


By neglecting the terms of truncation error from (3) and solving for w_i^{m+2} we get

(4)   \begin{equation*} w_i^{m+2}=4w_i^{m+1}-3w_i^m-\dfrac{2\alpha \Delta t}{\Delta x^2}[{w_{i+1}^m-2w_i^m+w_{i-1}^m}]\end{equation*}

m=0,1,\cdots,M-2,

where N=\dfrac{L}{\Delta x} and M=\dfrac{t_{max}}{\Delta t}.
Gerald W. Recktenwald denoted \dfrac{\alpha \Delta t}{\Delta x^2} by
\textit{r} ({r}=\dfrac{\alpha \Delta t}{\Delta x^2})
and is taken into
consideration for all of the schemes. Equation (4)
is known as
the FTCS approximation of the heat equation. After rearrangement we can write

(5)   \begin{equation*} w_i^{m+2} =4w_i^{m+1}-2r(w_{i-1}^m-2w_i^m+w_{i+1}^m)-3w_i^m \quad , m=0,1,\cdots,M-2, i=1,2,\cdots,N-1 \end{equation*}

BTCS

For deriving (4), the scheme of forward difference
has been applied
for approximating the time derivative at the left side of (??). Now, by choosing the scheme of backward difference,

(6)   \begin{equation*} \dfrac{\partial w}{\partial t}\big|{t{m},x_i}=\dfrac{3w_{i}^{m}-4w_{i}^{m-1}+w_{i}^{m-2}}{2\Delta t}+\mathcal{O}(\Delta t^2)\end{equation*}


Using the approximation of central difference for (\dfrac{\partial^2 w}{\partial x^2}){x_i} ,

(7)   \begin{equation*} \dfrac{\partial^2 w}{\partial x^2}\big|{x_i} =\dfrac{w_{i+1}^m-2w_i^m+w_{i-1}^m}{\Delta x^2}+\mathcal{O}(\Delta x^2)\end{equation*}


Using these equations (6) and (7) in (??) and combine the terms of truncation error to yield

(8)   \begin{equation*} \dfrac{3w_{i}^{m}-4w_{i}^{m-1}+w_{i}^{m-2}}{2\Delta t}=\alpha \big[\dfrac{w_{i+1}^m-2w_i^m+w_{i-1}^m}{\Delta x^2}\big]+\mathcal{O}(\Delta t^2 , \Delta x^2)\end{equation*}


Drop the truncation error terms from Equation (8) ,we can write

(9)   \begin{equation*} 3w_{i}^{m}-4w_{i}^{m-1}+w_{i}^{m-2}=\dfrac{2 \alpha \Delta t}{\Delta x^2} \big[{w_{i+1}^m-2w_i^m+w_{i-1}^m}\big] \end{equation*}


It can be written as

(10)   \begin{equation*} -\frac{\alpha }{\Delta x^2}w_{i-1}^m+(\frac{3}{2 \Delta t}+\frac{2 \alpha}{\Delta x^2})w_i^{m}-\frac{\alpha }{\Delta x^2}w_{i+1}^m=\frac{2}{\Delta t}w_i^{m-1}-\frac{1}{2 \Delta t}w_i^{m-2}\end{equation*}

\quad m=2,3,\cdots,M ,

Equation (10) is known as BTCS approximation for the heat equation.

Finite Difference Schemes for the Heat Eq: BTCS for two points
Finite Difference Schemes for the Heat Eq: BTCS for 2 and 5 points

As we know that Backward in time, centered in space is an implicit scheme and unconditionally stable.
So from the we can clearly see that, there is no inconsistency as shown in figures ?? and ??. It is clearly we can see that figure ?? the value of t=0.5 is maximum go to w(x,t)= 0.4 but in figure ?? this value of w attains the value about 0.7. So we can say that the scheme in backward in time, centered in space for three points in time discretizatoin has more refinement then two points.

Crank-Nicolson

In Crank-Nicolson method, we apply time backward difference already used for BTCS scheme to approximate the left side of heat equation. Whereas right side of the same equation have been approximated by
using average of central difference method calculated at the previous and current time step.
Hence, (??) is being approximated with,

(11)   \begin{equation*} \dfrac{3w_{i}^{m}-4w_{i}^{m-1}+w_{i}^{m-2}}{2\Delta t}=\dfrac{\alpha}{2} [\dfrac{w_{i+1}^m-2w_i^m+w_{i-1}^m}{\Delta x^2}+ \dfrac{w_{i+1}^{m-1}-2w_i^{m-1}+w_{i-1}^{m-1}}{\Delta x^2}]\end{equation*}


After some rearrangements we can write

(12)   \begin{equation*} 3w_i^m-r(w_{i-1}^m-2w_i^m+w_{i+1}^m)= r(w_{i-1}^{m-1}-2w_i^m+w_{i+1}^{m-1})+4w_i^{m-1}-w_i^{m-2}, \end{equation*}

m=2,3,4, \cdots ,M-1,M , \quad i=1,2,\cdots,N-2,N-1

Adams-Bashforth

From equation (??)

(13)   \begin{equation*} \dfrac{w_i^{m+1}-w_i^m}{\Delta t} = \dfrac{3g^m-g^{m-1}}{2} \end{equation*}


Since g=\dfrac{\partial w}{\partial t}= \alpha \dfrac{\partial^2 w}{\partial x^2}. So by using central difference for \dfrac{\partial^2 w}{\partial x^2}, we can write

(14)   \begin{equation*} w_i^{m+1}=w_i^m+\dfrac{r}{2}[3(w_{i-1}^m-2w_i^m+w_{i+1}^m)-(w_{i-1}^{m-1}-2w_i^{m-1}+w_{i+1}^{m-1})] \end{equation*}

Von Neumann Stability Analysis

Expanding Fourier series of a function w(x,t) in the spatial direction,

(15)   \begin{equation*} w(x, t) = \sum_{k=0}^{\infty} A_k(t) e^{ikx} \end{equation*}


where A_k(t) are the functions of time and k is the wave number.
To check the stability of a numerical scheme,
the amplitude A_k(t) must be bounded over time.
For the linear problems, instead of considering sum of all wave number components,
we can examine the stability of the amplitude for each
wave number k.
In terms of the amplitude A = \frac{A^{m+1}}{A^m}, the
conditions for the stability of numerical schemes are

(16)   \begin{align*} \big| A \big| \left\{ \begin{array}{cc} <1, & \hspace{5mm} \text{scheme is stable} \\ =1, & \hspace{3mm} \text{scheme is neutrally stable} \\ >1, & \hspace{5mm} \text{ scheme is unstable} \\ \end{array} \right. \end{align*}

Higher Order Finite Difference Schemes for the Heat Eq (Part 1)

Below we check the Von Neumann stability analysis of
some discretization schemes, discussed in the previous sections.

FTCS

For FTCS the Heat Equation can be discretized as in equation( ??)

(17)   \begin{equation*} w_i^{m+1}=w_i^m+ r (w_{i+1}^m-2w_i^m+w_{i-1}^m)\end{equation*}


where r=\frac{\alpha \Delta t}{ \Delta x^2}.
By using

(18)   \begin{equation*} w(x, t) = \sum_{k=0}^{\infty} A_k(t)e^{ikx} \end{equation*}


Equation (17) will become

(19)   \begin{equation*} A^{m+1}=A^m + r (e^{-ik \Delta x} - 2 + e^{ik \Delta x}) A^m\end{equation*}


Since cos(k \Delta x)= \frac{e^{-ik \Delta x} + e^{ik \Delta x}}{2} , so we can write

(20)   \begin{equation*} A^{m+1}=[1 + 2 r (cos (k \Delta x) -1)] A^m \end{equation*}


(21)   \begin{equation*} A^{m+1}=[1 - 4 r sin^2 (\dfrac{k \Delta x}{2})] A^m \end{equation*}


So, the amplitude A is

(22)   \begin{equation*} A=1 - 4 r sin^2 (\dfrac{k \Delta x}{2}) \end{equation*}


For this scheme to be stable , we require |A|\leq 1, which gives

    \[-1 \leq 1 - 4 r sin^2 (\dfrac{k \Delta x}{2})\leq 1\]


    \[4 r \leq 2\]


    \[\Rightarrow r \leq \dfrac{1}{2}\]

This analysis means that the FTCS scheme is conditionally stable
and the heat equation can be integrated stably with FTCS ,
by choosing the time step must be according to \Delta t \leq \dfrac{ \Delta x^2}{2 \alpha }.

Leave a Comment