# 1-D time dependent Parabolic differential equations

Parabolic differential equations could be written as:
$$\frac{\partial u(x,t)}{\partial t} + L[u(x,t)] = f(x,t), \ \ \ \ \ \ x \in [\alpha,\beta], \ t>0 \\$$
where $L$ is an elliptic operator applied to the solution, $f$ is a known function, and $u$ is the solution. Before implementing the finite difference schemes, it is important to consider some main features of the parabolic equations:
- The solution admits just one characteristic line;
- This line defines the dependence domain;
- Along the parabolic direction $t$ , just one “initial condition” is required;
- Each node at the edge needs a boundary condition;

A generic example of such a problem is:
$$\frac{\partial u(x,t)}{\partial t} = -c \frac{\partial u(x,t)}{\partial x} + \sigma \frac{\partial ^2 u(x,t)}{\partial x^2} + f(x,t);$$ and the more specific case is the diffusion equation:
$$\frac{\partial u(x,t)}{\partial t} = \sigma \frac{\partial ^2 u(x,t)}{\partial x^2} + f(x,t).$$ Because of the importance of this equation in describing physical phenomena, in the following it is considered as a paradigm of parabolic equations.

# Diffusion equation: finite difference discretization

Along the time direction, parabolicity of the equation suggests the use of forward scheme in order to consider the current solution as a function of just the previous time-steps.
Along the spatial direction, the homogeneous mathematical problem is symmetrical, so the best initial choice is to respect the symmetry; this way also guarantees the best convergence speed for constant $\Delta x$ grid.
Referring to the previous chapter, a first order forward scheme is used for the unsteady term, and the second order centred scheme is employed for the diffusion term as follows:
$$\frac{u^{n+1}-u^n}{\Delta t} = \sigma \frac{u_{i+1} - 2 u_i + u_{i-1}}{\Delta x^2} + f_i.$$

# Explicit Euler Scheme

This scheme is said to be the most convenient finite difference scheme for discretizing the ODEs. The diffusion equation can be discretized using this scheme as follows:
$$\frac{u_i^{n+1}-u_i^n}{\Delta t} = \sigma \frac{u_{i+1}^n - 2 u_i^n + u_{i-1}^n}{\Delta x^2} + f_i^n.$$
Fig. b.1. Explicit Euler scheme

Before writing the code for the discretization, it would be better to write down the discretized equation in matrix form:
$$\left[ I \right] \left(\{u\}^{n+1} - \{u\}^{n} \right) = \left[L\right] \{u\}^n + \Delta t \{f\}^n$$ where:
$$[L] = \left( \begin{array}{ccccc} -2\nu & \nu & 0 & \ldots & 0 \\ \nu & -2\nu & \nu & \ldots & 0 \\ 0 & \nu & -2\nu & \ldots & 0 \\ \vdots & & & \ddots & \vdots \\ 0 & 0 & 0 & \ldots & -2\nu \end{array} \right) \ , \ \nu = \frac{\sigma \Delta t}{\Delta x^2}$$ If $[S] = [I] + [L]$, then we have:
$$\{u\}^{n+1} = \left[S\right]\{u\}^{n} + \Delta t \{f\}^n$$

# Two Dirichlet boundary conditions

Imposing Dirichlet boundary conditions at each boundary of the diffusion equation means that a value of the function is specified on that boundary of the domain. Considering the matrix form of the diffusion equation this type of boundary condition is described in the following.
$$\left\lbrace\begin{array}{c} u_1\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1} \\ u_N \\ \end{array}\right\rbrace^{n+1} = \left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S] & & \\ & & & & \\ & & & & \\ & & & &\\ 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right) \left\lbrace\begin{array}{c} u_{Dl}\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1}\\ u_{Dr} \\ \end{array}\right\rbrace^n + \Delta t\left\lbrace\begin{array}{c} 0\\ f_2\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ 0 \\ \end{array}\right\rbrace^n$$ which can be also presented as:
$$\begin{cases} u^{n+1}_1 = u^{n+1}_{Dl} \\ u^{n+1}_i - u^{n}_i = \frac{\sigma \Delta t}{\Delta x^2} \left({u^n_{i-1} -2 u^n_i + u^n_{i+1}}\right) + f^n_i\Delta t \\ u^{n+1}_N = u^{n+1}_{Dr} \end{cases}$$ $u_{Dl}$ and $u_{Dr}$ are the value of Dirichlet boundary conditions at nodes $x_1$ and $x_N$ respectively. The highlighted parts of the vectors and the matrix show the implementation of the Dirichlet boundary conditions. For each time-step $n+1$, the values of $u$ are calculated using the given values of the previous time-step $n$.

# Neumann boundary condition on the left

In order to implement the Neumann boundary condition at node $x_1$, we have:
$$\left.\frac{\partial u}{\partial t}\right|_{x=x_1} = \frac{u_0^{n+1} - u_2^{n+1}}{2\Delta x} = u_{Nl};$$ where $u_{N1}$ is the value of Neumann boundary condition at node $x_1$. It means that a ghost node $x_0$ should be considered. However, using Eq. (5), $u_2^{n+1}$ can be obtained:
$$u_2^{n+1}-u_2^n = \nu \left(u_1^n - 2 u_2^n + u_3^n\right) + \Delta t f_2^n.$$ Consequently, we have:
$$u_0^{n+1} = \nu u_1^n + (1-2\nu) u_2^n + \nu u_3^n + 2 u_{Nl}^{n+1} \Delta x + f_2^n \Delta t.$$ As a result, the Neumann boundary condition at $x_1$ can be shown as following:
$$\left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_N \\ \end{array}\right\rbrace^{n+1} = \left( \begin{array}{cccccc} 0 & \nu & 1-2\nu & \nu & \ldots & 0 \\ & & & & & \\ & & & & & \\ & & [S] & & & \\ & & & & & \\ & & & & & \\ 0 & 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right) \left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_{Dr} \\ \end{array}\right\rbrace^n + \Delta t\left\lbrace\begin{array}{c} 2\frac{\Delta x}{\Delta t}u_{Nl} + f_2\\ f_1\\ \vdots\\ f_i\\ \vdots\\ 0 \\ \end{array}\right\rbrace^n$$

# Periodic boundary conditions

For the periodic boundary conditions, the last node is omitted from the solution vector. The following corrections on matrix $S$ describes the periodic boundary conditions:
$$[S] = \left( \begin{array}{cccccc} -2\nu+1 & \nu & 0 & \ldots & \nu & \not{0} \\ \nu & -2\nu+1 & \nu & \ldots & 0 & \not{0} \\ 0 & \nu & -2\nu+1 & \ldots & 0 & \not{0} \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ \nu & 0 & 0 & \ldots & -2\nu+1 & \not{\nu} \\ \not{0} & \not{0} & \not{0} & & \not{\nu} & \not{-}\not{2}\not{\nu}\not{+}\not{1}\\ \end{array} \right)$$ As it is shown above, the last column and row are eliminated.

# Implicit Euler scheme

This scheme is quite similar to the previous scheme, but differs in that it is an implicit method. The general form of this scheme is:
$$\frac{u^{n+1}_i - u^{n}_i}{\Delta t} = \sigma \frac{u^{n+1}_{i-1} -2 u^{n+1}_i + u^{n+1}_{i+1}}{\Delta x^2} + f^{n+1}_i$$
Fig. b.2. Implicit Euler scheme

The discretized equation in matrix form can be written as:
$$\left[ I \right] \left(\{u\}^{n+1} - \{u\}^{n} \right) = \left[L\right] \{u\}^{n+1} + \Delta t \{f\}^{n+1}$$ where:
$$[L] = \left( \begin{array}{ccccc} -2\nu & \nu & 0 & \ldots & 0 \\ \nu & -2\nu & \nu & \ldots & 0 \\ 0 & \nu & -2\nu & \ldots & 0 \\ \vdots & & & \ddots & \vdots \\ 0 & 0 & 0 & \ldots & -2\nu \end{array} \right) \ , \ \nu = \frac{\sigma \Delta t}{\Delta x^2}$$ If $[S] = [I] - [L]$, then we have:
$$\{u\}^{n+1} = \left[S\right]^{-1}\left(\{u\}^{n} + \Delta t \{f\}^{n+1}\right)$$

# Dirichlet boundary conditions

The Dirichlet boundary conditions on the left and right can be expressed as:
$$\left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S] & & \\ & & & & \\ & & & & \\ & & & &\\ 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right) \left\lbrace\begin{array}{c} u_1\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1} \\ u_N \\ \end{array}\right\rbrace^{n+1} = \left\lbrace\begin{array}{c} u_{Dl}\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1}\\ u_{Dr} \\ \end{array}\right\rbrace^n + \Delta t\left\lbrace\begin{array}{c} 0\\ f_2\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ 0 \\ \end{array}\right\rbrace^n$$ which can be also presented as:
$$\begin{cases} u^{n+1}_1 = u^{n+1}_{Dl} \\ u^{n+1}_i - u^{n}_i = \frac{\sigma \Delta t}{\Delta x^2} \left({u^{n+1}_{i-1} -2 u^{n+1}_i + u^{n+1}_{i+1}}\right) + f^{n+1}_i\Delta t \\ u^{n+1}_N = u^{n+1}_{Dr} \end{cases}$$

# Neumann boundary condition on the left

In order to implement the Neumann boundary condition at node $x_1$, we have:
$$\frac{u_0^{n+1} - u_2^{n+1}}{2\Delta x} = u^{n+1}_{Nl}$$ and in matrix form:
$$\left( \begin{array}{ccc} 0 & \ldots & 0 \\ & & \\ & & \\ [I] \\ & & \\ & &\\ & & \\ \end{array} \right) \left(\left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_N \\ \end{array}\right\rbrace^{n+1} - \left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_N \\ \end{array}\right\rbrace^{n} \right) = \left( \begin{array}{c} -1 \ \ 0 \ \ 1 \ \ \ldots \ \ 0 \\ \\ \\ [L] \\ \\ \\ \\ \end{array} \right) \left\lbrace\begin{array}{c} u_{0}\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_{N} \\ \end{array}\right\rbrace^{n+1} + \Delta t\left\lbrace\begin{array}{c} 2\frac{\Delta x}{\Delta t}u_{Nl}\\ f_1\\ \vdots\\ f_i\\ \vdots\\ f_{N} \\ \end{array}\right\rbrace^{n+1}$$ It should be mentioned that in this case, the right boundary condition is not set.

# Periodic boundary conditions

For the periodic boundary conditions, the last node is omitted from the solution vector. The following corrections on matrix $S$ describes the periodic boundary conditions:
$$[S] = \left( \begin{array}{cccccc} 2\nu+1 & -\nu & 0 & \ldots & -\nu & \not{0} \\ -\nu & 2\nu+1 & -\nu & \ldots & 0 & \not{0} \\ 0 & -\nu & 2\nu+1 & \ldots & 0 & \not{0} \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ -\nu & 0 & 0 & \ldots & 2\nu+1 & \not{-}\not{\nu} \\ \not{0} & \not{0} & \not{0} & & \not{-}\not{\nu} & \not{2}\not{\nu}\not{+}\not{1}\\ \end{array} \right)$$

# Crank-Nicolson scheme

The Crank-Nicolson scheme can be presented as follows:
$$\frac{u^{n+1}_i - u^{n}_i}{\Delta t} = \frac{1}{2}\sigma\left( \frac{u^n_{i-1} -2 u^n_i + u^n_{i+1}}{\Delta x^2} + f^n_i\right) + \frac{1}{2}\sigma \left(\frac{u^{n+1}_{i-1} -2 u^{n+1}_i + u^{n+1}_{i+1}}{\Delta x^2} + f^{n+1}_i\right)$$
Fig. b.3. Crank-Nicolson scheme

In order to write the method in matrix form, it is better to divide the explicit ($EE$) and implicit ($IE$) terms:
$$u^{n+1}_i - u^{n}_i = \frac{1}{2} rhs_{EE} + \frac{1}{2} rhs_{IE}$$ So we have:
$$\left[ I \right] \left(\{u\}^{n+1} - \{u\}^{n} \right) = \frac{1}{2}\left[L_{EE}\right] \{u\}^{n} + \frac{1}{2}\left[L_{IE}\right] \{u\}^{n+1} + \Delta t \left(\frac{1}{2}\{f\}^{n}+\frac{1}{2}\{f\}^{n+1}\right)$$ If $[S_{EE}] = [I] + \frac{1}{2}[L_{EE}]$ and $[S_{IE}] = [I] - \frac{1}{2}[L_{IE}]$ we can rewrite the Eq. 27 as:
$$\{u\}^{n+1} = \left[S_{IE}\right]^{-1} \left(\left[S_{EE}\right]\{u\}^{n} + \frac{\Delta t}{2} \left(\{f\}^{n}+\{f\}^{n+1}\right)\right)$$

# Dirichlet boundary conditions

$$\left( \begin{array}{ccccc} 1 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S_{IE}] & & \\ & & & & \\ & & & & \\ & & & &\\ 0 & 0 & 0 & \ldots & 1 \\ \end{array} \right) \left\lbrace\begin{array}{c} u_1\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1} \\ u_N \\ \end{array}\right\rbrace^{n+1} = \left( \begin{array}{ccccc} 0 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S_{EE}] & & \\ & & & & \\ & & & & \\ & & & &\\ 0 & 0 & 0 & \ldots & 0 \\ \end{array} \right)\left\lbrace\begin{array}{c} u_1\\ u_2\\ \vdots\\ u_i\\ \vdots\\ u_{N-1}\\ u_N \\ \end{array}\right\rbrace^n + \frac{\Delta t}{2}\left(\left\lbrace\begin{array}{c} \frac{2 u_{Dl}}{\Delta t}\\ f_2\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ \frac{2 u_{Dr}}{\Delta t} \\ \end{array} \right\rbrace^n + \left\lbrace\begin{array}{c} 0\\ f_2\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ 0 \\ \end{array}\right\rbrace^{n+1}\right)$$

# Neumann boundary condition on the left

$$\left( \begin{array}{ccccc} 1 & 0 & -1 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S_{IE}] & & \\ & & & & \\ & & & & \\ & & & &\\ & & & & \\ \end{array} \right) \left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_{N-1} \\ u_N \\ \end{array}\right\rbrace^{n+1} = \left( \begin{array}{ccccc} 0 & 0 & 0 & \ldots & 0 \\ & & & & \\ & & & & \\ & & [S_{EE}] & & \\ & & & & \\ & & & & \\ & & & &\\ & & & & \\ \end{array} \right)\left\lbrace\begin{array}{c} u_0\\ u_1\\ \vdots\\ u_i\\ \vdots\\ u_{N-1}\\ u_N \\ \end{array}\right\rbrace^n + \frac{\Delta t}{2}\left(\left\lbrace\begin{array}{c} \frac{4\Delta x u_{Nl}}{\Delta t}\\ f_1\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ f_{N} \\ \end{array} \right\rbrace^n + \left\lbrace\begin{array}{c} 0\\ f_1\\ \vdots\\ f_i\\ \vdots\\ f_{N-1}\\ f_{N} \\ \end{array}\right\rbrace^{n+1}\right)$$

# Periodic boundary conditions

$$[S_{IE}] = \left( \begin{array}{cccccc} 2\nu+1 & -\nu & 0 & \ldots & -\nu & \not{0} \\ -\nu & 2\nu+1 & -\nu & \ldots & 0 & \not{0} \\ 0 & -\nu & 2\nu+1 & \ldots & 0 & \not{0} \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ -\nu & 0 & 0 & \ldots & 2\nu+1 & \not{-}\not{\nu} \\ \not{0} & \not{0} & \not{0} & & \not{-}\not{\nu} & \not{2}\not{\nu}\not{+}\not{1}\\ \end{array} \right)$$
$$[S_{EE}] = \left( \begin{array}{cccccc} -2\nu+1 & \nu & 0 & \ldots & \nu & \not{0} \\ \nu & -2\nu+1 & \nu & \ldots & 0 & \not{0} \\ 0 & \nu & -2\nu+1 & \ldots & 0 & \not{0} \\ & & & & & \\ & & & & & \\ & & & & & \\ & & & & & \\ \nu & 0 & 0 & \ldots & -2\nu+1 & \not{\nu} \\ \not{0} & \not{0} & \not{0} & & \not{\nu} & \not{-}\not{2}\not{\nu}\not{+}\not{1}\\ \end{array} \right)$$

The relevant Matlab codes can be downloaded using the following link:
Matlab Codes Bank