Skip to main content

Section 2.4 Finite Differences Method: Parabolic PDEs

Just as the advection PDE, also the heat PDE comes from the continuity equation
\begin{equation*} \frac{\partial u}{\partial t} + {\mathbf\nabla}\cdot\mathbf j = 0. \end{equation*}
The difference is that here we do not consider a fluid, where \(\mathbf j=u{\mathbf v}\text{,}\) but rather some diffusing material that follows Fick's law, namely

\begin{equation*} \mathbf j(t,{\mathbf x}) = -K(u(t,{\mathbf x}),{\mathbf x}){\mathbf\nabla} u(t,{\mathbf x}), \end{equation*}

where the function \(K\) is the diffusion coefficient.

The resulting PDE is called Diffusion Equation:
\begin{equation*} u_t = \mathbf\nabla\cdot (K\mathbf \nabla u). \end{equation*}

Heat equation. When the diffusion coefficient is constant, the equation is called Heat PDE
\begin{equation*} u_t = K\Delta u, \end{equation*}
where \(\Delta\) is the Laplacian operator, namely
\begin{equation*} \Delta u(x^1,\dots,x^n) = u_{x^1x^1}+\dots+u_{x^nx^n}. \end{equation*}
1D Heat equation. In 1+1 dimensions, the heat PDE writes
\begin{equation*} u_t=c\,u_{xx}\,,\;u(0,x)=f(x)\,. \end{equation*}
We consider Dirichlet boundary conditions \(u(t,0)=u(t,1)=0.\) Physically, this corresponds to heating a wire with some initial temperature distribution \(f(x)\) while keeping the end-points at some fixed temperature.

The solution of the PDE will give the evolution of the temperature as \(t\to\infty\text{.}\) Solutions to
\begin{equation*} u_t=c\,u_{xx},\;u(0,x)=f(x),u(t,0)=u(t,1)=0, \end{equation*}
with separate variables can be found by setting \(u(t,x)=T(t)X(x)\text{,}\) leading to the ODE pair
\begin{equation*} T'=\lambda cT,\;X'=\lambda X^2. \end{equation*}
By enforcing the boundary conditions, we see that \(X(0)=X(1)=0\) and so \(\lambda=-(k\pi)^2\text{.}\)

In general, every solution can be written in Fourier series as
\begin{equation*} u(t,x)=\sum_{k=1}^\infty A_k e^{-c(k\pi)^2t}\sin(k\pi x), \end{equation*}
with
\begin{equation*} A_k=2\int_0^1f(\xi)\sin(k\pi\xi)d\xi. \end{equation*}

Subsection 2.4.1 FTCS

As in case of the advection, we use the approximations
\begin{equation*} u_t(t_n,x_j)\to\frac{u_{n+1,j}-u_{n,j}}{\Delta t} \end{equation*}
and
\begin{equation*} u_{xx}(t_n,x_j)\to\frac{u_{n,j+1}-2u_{n,j}+u_{n,j-1}}{\left(\Delta x\right)^2}, \end{equation*}
so the equation translates into
\begin{equation*} u_{n+1,j} = u_{n,j} + \nu\left(u_{n,j+1}-2u_{n,j}-u_{n,j-1}\right), \end{equation*}
where \(\nu=\frac{c\Delta t}{\left(\Delta x\right)^2}\)
FTCS Stability. As for the hyperbolic case, we consider the solution \(u(t,x) = e^{ik(x-c t)}\text{,}\) so that
\begin{equation*} u_{n,j} = e^{ik(j\Delta x-c n\Delta t)} = G(k)^n e^{ikj\Delta x}\,, \end{equation*}
and we plug it into the system, so that we get
\begin{equation*} G(k)^{n+1} e^{ikj\Delta x} = G(k)^{n} e^{ikj\Delta x} + \nu\left(G(k)^{n} e^{ik(j+1)\Delta x} - 2G(k)^{n} e^{ikj\Delta x} + G(k)^{n} e^{ik(j-1)\Delta x}\right)\,, \end{equation*}
namely
\begin{equation*} G(k) = 1 - \nu(2\cos\left(k\Delta x\right)-2) = 1 - 4\nu\sin^2\frac{k\Delta x}{2}. \end{equation*}
Hence, necessary condition for the stability of this method is
\begin{equation*} \nu=\frac{c\Delta t}{(\Delta x)^2}\leq\frac{1}{2} \end{equation*}

Remark 2.4.1.
The FTCS method has a big problem: in order to keep it stable, if we decrease \(\Delta x\) by a factor \(h\) we need to decrease at the same time \(\Delta t\) by a factor \(h^2\text{,}\) which leads to unreasonably small values of \(\Delta t\) when a large accuracy is needed.

Subsection 2.4.2 CTCS (Richardson Method)

One can try using central derivative approx for time:
\begin{equation*} u_t(t_n,x_j)\to\frac{u_{n+1,j}-u_{n-1,j}}{2\Delta t}. \end{equation*}
Then the system becomes
\begin{equation*} u_{n+1,j} = u_{n-1,j} + 2\nu\left(u_{n,j+1}-2u_{n,j}-u_{n,j-1}\right)\,. \end{equation*}
The von Neumann analysis shows that the amplification factor \(G(k)\) satisfies
\begin{equation*} g^2+2\beta g-1,\;\;\beta=4\nu\sin^2\frac{k\Delta x}{2}, \end{equation*}
which has a solution with modulus larger than 1 for all \(\nu\text{.}\) Hence, this method is unconditionally unstable.

Subsection 2.4.3 DuFort-Frankel method

Similarly to how we did for the Lax algorithm of the Advection PDE, one can try replacing the \(u_{n,j}\) term with its time average:
\begin{equation*} u_{n+1,j} = u_{n-1,j} + 2\nu\left(u_{n,j+1}-(u_{n+1,j}+u_{n-1,j})+u_{n,j-1}\right). \end{equation*}
Hence
\begin{equation*} u_{n+1,j}=\frac{1-2\nu}{1+2\nu}u_{n-1,j}+2\frac{\nu}{1+2\nu}\left(u_{n,j+1}+u_{n,j-1}\right). \end{equation*}
Now \(G(k)\) solves
\begin{equation*} (1+2\nu)g^2+4g\nu\cos(k\Delta x)+2\nu-1=0\, \end{equation*}
whose roots are all outside of the unit disk for all \(\nu\text{,}\) namely this method is unconditionally stable!

Does this method converge to the exact solution? Actually not! Indeed, while it is stable, it is not consistent.

We can verify this with Taylor expansion. Write just \(u_t\text{,}\) \(u_x\) and so on for the values of those functions at \((t_n,x_j)\text{.}\) Then:
\begin{equation*} u_{n+1,j}-u_{n,j} = u_t \Delta t + O(\Delta t^3), \end{equation*}
\begin{equation*} u_{n,j+1}-(u_{n+1,j}+u_{n-1,j})+u_{n,j-1} = u_{xx}\Delta x^2 - u_{tt}\Delta t^2 + O(\Delta x^4,\Delta t^4) \end{equation*}
and so the DuFort-Frankel method system is equivalent to
\begin{equation*} u_t = c u_{xx} - c \frac{\Delta t^2}{\Delta x^2} u_{tt} + O(\Delta x^2,\Delta t^2,\Delta t^4/\Delta x^2). \end{equation*}
This shows that this method is not consistent for all choices of \(\Delta t,\Delta x\to0\text{.}\)

For instance, if we choose \(\Delta t=\Delta x=h\text{,}\) then the numerical solutions of the DuFort-Frankel method actually converge, for \(h\to0\text{,}\) to the solutions of the hyperbolic PDE
\begin{equation*} u_t = c ( u_{xx} - u_{tt} )\;\;! \end{equation*}

Subsection 2.4.4 BTCS method

Often implicit methods fix many problems, so one can try replacing the time derivative with its backward in time expression:
\begin{equation*} u_t(t_n,x_j)\to\frac{u_{n,j}-u_{n-1,j}}{\Delta t} \end{equation*}
The system now becomes

\begin{equation*} u_{n,j} = u_{n-1,j} + \nu\left(u_{n,j+1}-2u_{n,j}+u_{n,j-1}\right)\,, \end{equation*}

that can be solved "backwards":

\begin{equation*} u_{n-1,j} = u_{n,j} - \nu\left(u_{n,j+1}-2u_{n,j}+u_{n,j-1}\right)\,, \end{equation*}
Namely, \(\vec u_{n-1}\text{,}\) the vector of values of \(u\) at \(t=t_{n-1}\text{,}\) is equal to \(M\vec u_{n}\text{,}\) where the matrix \(M\) is tridiagonal with \(1+2\nu\) on the diagonal and \(-\nu\) in the upper/lower diagonal terms.

Hence
\begin{equation*} \vec u_{n} = (M^{-1})^n\cdot\vec u_{0}. \end{equation*}

Stability. For BTCS we get
\begin{equation*} G(k)^{n+1} e^{ikj\Delta x} = G(k)^{n} e^{ikj\Delta x} + \nu\left(G(k)^{n+1} e^{ik(j+1)\Delta x} - 2G(k)^{n+1} e^{ikj\Delta x} + G(k)^{n+1} e^{ik(j-1)\Delta x}\right), \end{equation*}
namely
\begin{equation*} G(k) = 1 - \nu\,G(k)(2\cos\left(k\Delta x\right)-2), \end{equation*}
so that
\begin{equation*} G(k) = \frac{1}{1+4\nu\sin^2\frac{k\Delta x}{2}}. \end{equation*}
This method therefore is stable for all \(\Delta t\) and \(\Delta x\text{!}\)

A weak point. Via Taylor series analysis, we see that the BTCS system is equivalent to
\begin{equation*} u_t = c u_{xx} + O(\Delta t,\Delta x^2). \end{equation*}
Just as with FTCS, being only order 1 in time is a problem when a reasonably small error in numerical solutions is needed.

Implementation in MATLAB.

Subsection 2.4.5 Crank-Nicholson method

Replace the CS r.h.s. with its forward average in time:
\begin{equation*} u_{xx}(t_n,x_j)\to \frac{u_{n+1,j+1}-2u_{n+1,j}+u_{n+1,j-1}}{\left(\Delta x\right)^2} + \frac{u_{n,j+1}-2u_{n,j}+u_{n,j-1}}{\left(\Delta x\right)^2} \end{equation*}
The Heat equation now translates into
\begin{equation*} -\nu u_{n+1,j+1}+2(1+\nu)u_{n+1,j}-\nu u_{n+1,j-1} = \nu u_{n,j+1}+2(1-\nu)u_{n,j}+\nu u_{n,j-1} \end{equation*}
which writes as \(M\mathbf u_{n+1}=N\mathbf u_{n}\) for two tridiagonal matrices \(M,N\text{,}\) so that
\begin{equation*} \mathbf u_{n} = (M^{-1}N)^n\cdot\mathbf u_{0} \end{equation*}
In this case we get
\begin{equation*} G(k) = \frac{1-2r\sin^2\frac{k\Delta x}{2}}{1+2r\sin^2\frac{k\Delta x}{2}} \end{equation*}

This method therefore is also both unconditionally stable and second order in both time and space.

Implementation in MATLAB

Subsection 2.4.6 The 2D Heat PDE

Below we solve the Heat PDE
\begin{equation*} u_{t} = \alpha\left(u_{xx}+u_{yy}\right). \end{equation*}
Extending the method above to the case of 2 space dimensions, we find the algorithm
\begin{equation*} u_{n+1,m,l} = u_{n,m,l} + \alpha\frac{\Delta t}{\Delta x^2}\left(u_{n,m+1,l}+u_{n,m,l+1}-4u_{n,m,l}+u_{n,m-1,l}+u_{n,m,l-1}\right). \end{equation*}
The coefficient of \(u_{n,m,l}\) must be positive to have stability, so we must ask
\begin{equation*} \frac{\Delta x^2}{\Delta t}\geq4\alpha \end{equation*}
The code below solves the wave equation above in the square \(S=[0,1]\times[0,1]\) with initial conditions
\begin{equation*} u(0,x,y) = e^{-100((x-0.5)^2+(y-0.5)^2)},\;u_t(0,x,y)=0 \end{equation*}
and boundary conditions \(u_{\partial S}=0\text{.}\)