function to compute the right hand side of heat equation
function rhs(nx,dx,u,r,alpha)
for i = 2:nx
r[i] = alpha*(u[i+1] - 2.0*u[i] + u[i-1])/(dx*dx)
end
end
```
✝ ✆

**Figure 3.** Comparison of exact solution and numerical solution computed using third-order Runge-Kutta numerical scheme. The numerical solution is computed using Δ*t* = 0.0025 and Δ*x* = 0.025.

#### *2.3. Crank-Nicolson Scheme*

The Crank-Nicolson scheme is a second-order accurate in time numerical scheme. The Crank-Nicolson scheme is derived by combining forward in time method at (*n*) and the backward in time method at (*n* + 1). Similar to Equation (3), we can write the Taylor series in time for the discrete field *<sup>u</sup>*(*n*+1) *<sup>i</sup>* as

$$u\_i^{(n)} = u\_i^{(n+1)} - \frac{\partial u\_i^{(n+1)}}{\partial t} \Delta t + \mathcal{O}(\Delta t^2),\tag{15}$$

$$\frac{\partial u\_i^{(n+1)}}{\partial t} = \frac{u\_i^{(n+1)} - u\_i^{(n)}}{\Delta t} + \mathcal{O}(\Delta t). \tag{16}$$

The second-derivative term in heat equation at (*n* + 1)*th* time step can be approximated using the second-order central difference scheme similar to the way we derived in Section 2.1

$$\frac{\partial^2 u\_i^{(n+1)}}{\partial \mathbf{x}^2} = \frac{u\_{i+1}^{(n+1)} - 2u\_i^{(n+1)} + u\_{i-1}^{(n+1)}}{\Delta \mathbf{x}^2} + \mathcal{O}(\Delta \mathbf{x}^2). \tag{17}$$

Substituting finite difference approximation derived in Equations (16) and (17) in the heat equation, we get

$$\frac{u\_i^{(n+1)} - u\_i^{(n)}}{\Delta t} = a \frac{u\_{i+1}^{(n+1)} - 2u\_i^{(n+1)} + u\_{i-1}^{(n+1)}}{\Delta x^2}.\tag{18}$$

If we add Equations (8) and (18), the leading term in discretization error cancels each other due to their opposite sign and same magnitude. Hence, we get the second-order accurate in time numerical scheme. The Crank-Nicolson scheme is O(Δ*t* 2, Δ*x*2) numerical scheme. The Crank-Nicolson scheme is given as

$$\frac{u\_i^{(n+1)} - u\_i^{(n)}}{\Delta t} = \frac{u}{2} \left[ \frac{u\_{i+1}^{(n)} - 2u\_i^{(n)} + u\_{i-1}^{(n)}}{\Delta x^2} + \frac{u\_{i+1}^{(n+1)} - 2u\_i^{(n+1)} + u\_{i-1}^{(n+1)}}{\Delta x^2} \right]. \tag{19}$$

We observe from the above equation that the solution at (*n* +1)*th* time step depends on the solution at previous time step (*n*)*th* and the solution at (*n* + 1)*th* time step. Such numerical schemes are called implicit numerical schemes. For implicit numerical schemes, the system of algebraic equations has to be solved by inverting the matrix. For one-dimensional problems, the tridiagonal matrix is formed and can be solved efficiently using the Thomas algorithm, which gives a fast O(*N*) direct solution as opposed to the usual <sup>O</sup>(*N*3) for a full matrix, where *<sup>N</sup>* is the size of the square tridiagonal matrix. The main advantage of an implicit numerical scheme is that they are unconditionally stable. This allows the use of larger time step and grid size without affecting the convergence of the numerical scheme. An illustration of a tridiagonal system is given in Equation (20). The Thomas algorithm [17] used for solving the tridiagonal matrix is outlined in Algorithm 1.

$$
\begin{bmatrix} b\_1 & c\_1 & 0 & \cdots & 0 \\ a\_2 & b\_2 & c\_2 & \cdots & 0 \\ \vdots & \vdots & & \ddots & \vdots \\ 0 & \cdots & a\_{N-1} & b\_{N-1} & c\_{N-1} \\ 0 & 0 & \cdots & a\_N & b\_N \end{bmatrix} \begin{bmatrix} u\_1 \\ u\_2 \\ \vdots \\ u\_{N-1} \\ u\_N \end{bmatrix} = \begin{bmatrix} d\_1 \\ d\_2 \\ \vdots \\ d\_{N-1} \\ d\_N \end{bmatrix} \tag{20}
$$



The implementation of the Crank-Nicolson numerical scheme for the heat equation is presented in Listing 3 and Thomas algorithm in Listing 4. The comparison of exact and numerical solution computed using the Crank-Nicolson scheme, and absolute error plot are displayed in Figure 4. The absolute error between the exact and numerical solution is less for the Crank-Nicolson scheme than the FTCS numerical scheme due to smaller truncation error.

**Listing 3.** Implementation of Crank-Nicolson numerical scheme for heat equation in Julia. ✞ ☎

✝ ✆

```