s,e: start and end index of the matrix
function tdms(a,b,c,r,x,s,e)
gam = Array{Float64}(undef, e)
bet = b[s]
x[s] = r[s]/bet
for i = s+1:e
gam[i] = c[i-1]/bet
bet = b[i] - a[i]*gam[i]
x[i] = (r[i] - a[i]*x[i-1])/bet
end
for i = e-1:-1:s
x[i] = x[i] - gam[i+1]*x[i+1]
end
return x
end
✝ ✆
```
#### *2.4. Implicit Compact Pade (ICP) Scheme*

We usually use the order of accuracy as an indicator of the ability of finite difference schemes to approximate the exact solution as it tells us how the discretization error will reduce with mesh refinement. Another way of measuring the order of accuracy is the modified wavenumber approach. In this approach, we see how much the modified wave number is different from the true wave number. The solution of a nonlinear partial differential equation usually contains several frequencies and the modified wavenumber approach provides a way to assess how well the different components of the solution are represented. The modified wavenumber varies for every finite difference scheme and can be found using Fourier analysis of the differencing scheme [12].

Compact finite difference schemes have very good resolution characteristics and can be used for capturing high-frequency waves [18]. In compact formulation, we express the derivative of a function as a linear combination of values of function defined on a set of nodes. The compact method tries to mimic global dependence and hence has good resolution characteristics. Lele [18] investigated the behavior of compact finite difference schemes for approximating a first and second derivative. The fourth-order accurate approximation to the second derivative is given by

$$\frac{1}{12}\frac{\partial^2 u}{\partial \mathbf{x}^2}\Big|\_{i-1} + \frac{10}{12}\frac{\partial^2 u}{\partial \mathbf{x}^2}\Big|\_i + \frac{1}{12}\frac{\partial^2 u}{\partial \mathbf{x}^2}\Big|\_{i+1} = \frac{u\_{i-1} - 2u\_i + u\_{i+1}}{\Delta \mathbf{x}^2}.\tag{21}$$

Taking the linear combination of heat equation and using the same coefficients as the above equation, we get

$$\frac{1}{12}\frac{\partial u}{\partial t}\Big|\_{i-1} + \frac{10}{12}\frac{\partial u}{\partial t}\Big|\_{i} + \frac{1}{12}\frac{\partial u}{\partial t}\Big|\_{i+1} = \frac{u\_{i-1} - 2u\_i + u\_{i+1}}{\Delta x^2}.\tag{22}$$

We then use Crank-Nicolson numerical scheme for time discretization and we arrive to below equation

$$\frac{u\_{i-1}^{(n+1)} - u\_{i-1}^{(n)}}{12\Delta t} + 10\frac{u\_i^{(n+1)} - u\_i^{(n)}}{12\Delta t} + \frac{u\_{i+1}^{(n+1)} - u\_{i+1}^{(n)}}{12\Delta t} = \frac{u\_{i-1}^{(n+1)} - 2u\_i^{(n+1)} + u\_{i+1}^{(n+1)} + u\_{i-1}^{(n)} - 2u\_i^{(n)} + u\_{i+1}^{(n)}}{2\Delta x^2} \tag{23}$$

Simplifying above equation, we get the implicit compact Pade (ICP) scheme for heat equation

$$a\_i u\_{i-1}^{(n+1)} + b\_i u\_i^{(n+1)} + c\_i u\_{i+1}^{(n+1)} = r\_i^{(n)},\tag{24}$$

where

$$a\_i = \frac{12}{\Delta x^2} - \frac{2}{a\Delta t'}\tag{25}$$

$$b\_i = \frac{-24}{\Delta x^2} - \frac{20}{a\Delta t},\tag{26}$$

$$c\_i = \frac{12}{\Delta x^2} - \frac{2}{a\Delta t'}\tag{27}$$

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

The ICP numerical scheme forms the tridiagonal matrix which can be solved using the Thomas algorithm to get the solution. The implementation of the ICP numerical scheme in Julia is given in Listing 5. Figure 5 shows the comparison of the exact and numerical solution and discretization error for ICP scheme. The ICP scheme is O(Δ*t* 2, Δ*x*4) accurate scheme and hence, the discretization error is very small.

**Listing 5.** Implementation of implicit compact scheme Pade numerical scheme for heat equation in Julia. ✞ ☎

```