Manoj Solanki

Project Abstract

A higher order numerical approach for solving 1-D parabolic PDEs, such as the Diffusion and Advection-Diffusion-Reaction (ADR) equations, is developed and studied in this work. The proposed approach uses theCrank–Nicolson traditional finite difference method (FDM) for time discretization and a fourth-order compactfinite difference method (CFDM) for spatial discretization. The accuracy and efficiency limits of the explicit,implicit, and Crank-Nicolson finite difference methods (FDM) are highlighted in this we use operator-based and Taylor series techniques, a compact study is created to overcome these issues, producing a solution that uses fewer grid points while achieving higher accuracy. The stability of the scheme is analyzed using Von Neumann ( Fourier analysis). A comparative study of CFDM is conducted through Python implementation on a variety of problems involving diffusion, advection, and reaction with constant coefficients. The results show that CFDM yields more accurate and efficient solutions than traditional methods. The work concludes by outlining future directions, particularly the extension of the scheme to PDEs with variable coefficients, for which initial developments are under way.

Introduction

Partial Differential Equations (PDEs) arise naturally in the mathematical modeling of physical, biological, and engineering phenomena involving spatial and temporal variations. Processes such as heat conduction, mass diffusion, fluid transport, chemical reactions, and population dynamics are commonly governed by parabolic PDEs. Due to the complexity of real-world systems, analytical solutions are often unavailable, making numerical methods an essential tool for practical simulations.

Among numerical techniques, finite difference methods (FDMs) are particularly attractive because of their conceptual simplicity, structured formulation, and ease of implementation on uniform grids. Classical methods such as the Forward-Time Central-Space (FTCS), Backward-Time Central-Space (BTCS), and Crank–Nicolson schemes are widely used for solving parabolic problems. However, these methods are typically limited to second-order accuracy in space, which may be insufficient for problems involving sharp gradients or advection-dominated behavior.

To address these limitations, compact finite difference methods (CFDMs) have been developed. These methods achieve higher-order spatial accuracy while using compact stencils, thereby improving computational efficiency without increasing stencil width or significantly adding to computational complexity. As a result, CFDMs provide a powerful and accurate alternative for the numerical solution of parabolic partial differential equations encountered in science and engineering applications.

Parabolic Partial Differential Equations

This thesis focuses on the numerical solution of the one-dimensional parabolic equation of the form

\[ \frac{\partial u}{\partial t} +\beta\frac{\partial u}{\partial x} +\gamma u = \alpha\frac{\partial^2u}{\partial x^2}, \qquad x\in(0,L),\; t>0, \tag{1} \]

where

u(x,t) represents the unknown scalar field,
$\alpha > 0$ is the diffusion coefficient,
$\beta$ is the advection velocity, and
$\gamma$ is the reaction coefficient.

Depending on the values of $\beta$ and $\gamma$, this equation reduces to:
The diffusion equation ($\beta = 0,\gamma = 0$),
The advection–diffusion equation ($\beta \neq 0,\gamma = 0$),
The advection–diffusion–reaction equation ($\beta \neq 0,\gamma \neq 0$),
These models arise in numerous applications such as heat transfer, pollutant transport, chemical kinetics, and biological systems.

Motivation for Compact Finite Difference Methods

Traditional finite difference schemes approximate spatial derivatives using Taylor series expansions truncated at second order. While simple, this approach requires extremely fine grids to achieve high accuracy, leading to increased computational cost and memory usage. Compact finite difference methods overcome this drawback by implicitly incorporating higher-order derivative information through algebraic manipulation of the governing equation itself. As a result:

fourth-order spatial accuracy is achieved,
only nearest-neighbor grid points are used,
numerical dispersion and artificial diffusion are significantly reduced,
This makes CFDMs particularly suitable for advection-dominated and stiff problems, where standard methods may produce oscillatory or overly diffusive solutions.

Main Contributions of This Thesis

The principal contributions of this dissertation are:

Development of an operator-based fourth-order compact scheme for the diffusion equation.
Construction of a general fourth-order compact finite difference scheme applicable to ADR equations.
Combination of compact spatial discretization with the Crank–Nicolson time integration.
Rigorous Von Neumann stability analysis of the proposed schemes.
Extensive numerical validation demonstrating fourth-order spatial convergence.
Python-based implementation and comparative analysis with classical schemes.

Main Scheme and Derivation

$$ Fu_{i+1}^{j+1}+Gu_i^{j+1}+Hu_{i-1}^{j+1} = \tilde{F}u_{i+1}^{j} +\tilde{G}u_i^{j} +\tilde{H}u_{i-1}^{j} \tag{2} $$

Where,

$$ F=-\frac{(K_1+K_2)}{2}-\frac{(K_3+K_4)}{k}, \qquad \tilde{F}=\frac{(K_1+K_2)}{2}-\frac{(K_3+K_4)}{k} $$ $$ G=\frac{2K_1-\gamma}{2}-\frac{(1-2K_3)}{k}, \qquad \tilde{G}=-\frac{2K_1-\gamma}{2}-\frac{(1-2K_3)}{k} $$ $$ H=-\frac{(K_1-K_2)}{2}-\frac{(K_3-K_4)}{k}, \qquad \tilde{H}=\frac{(K_1-K_2)}{2}-\frac{(K_3-K_4)}{k} $$ $$ K_1= -\left( \frac{\alpha}{h^2} +\frac{\beta^2}{12\alpha} -\frac{\gamma}{12} \right), \qquad K_2= \frac{\beta}{2h} -\frac{\gamma\beta h}{24\alpha}, \qquad K_3=\frac{1}{12}, \qquad K_4=-\frac{\beta h}{24\alpha} $$

For $i=1,2,\ldots,m-1$ and $j=0,1,\ldots,n$ we have,

$$ Fu_2^{j+1}+Gu_1^{j+1}+Hu_0^{j+1} = \tilde{F}u_2^j+\tilde{G}u_1^j+\tilde{H}u_0^j $$ $$ \vdots $$ $$ Fu_m^{j+1}+Gu_{m-1}^{j+1}+Hu_{m-2}^{j+1} = \tilde{F}u_m^j+\tilde{G}u_{m-1}^j+\tilde{H}u_{m-2}^j \tag{3} $$

with the help of the boundary condition, we get the following two relations

$$ u^{j+1}(x_0)=\Psi_L(t^{j+1}) \tag{4} $$ $$ u^{j+1}(x_m)=\Psi_R(t^{j+1}) \tag{5} $$

The resulting equation (3) for $i=2,3,\ldots,m-1$ yields a system of $m-1$ unknowns $u_1^{j+1},u_2^{j+1},\ldots,u_{m-1}^{j+1}$ at the $(j+1)^{th}$ time level for the fully discrete scheme for ADR equation which can be written in matrix form as follows,

$$ AU^{j+1}=BU^j+C \tag{6} $$

where,

$$ A= \begin{bmatrix} G & F & 0 & 0 & \cdots & 0 & 0 \\ H & G & F & 0 & \cdots & 0 & 0 \\ 0 & H & G & F & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & H & G & F & 0 \\ 0 & 0 & \cdots & 0 & H & G & F \\ 0 & 0 & \cdots & 0 & 0 & H & G \end{bmatrix}_{(m-1)\times(m-1)} $$ $$ B= \begin{bmatrix} \tilde{G} & \tilde{F} & 0 & 0 & \cdots & 0 & 0 \\ \tilde{H} & \tilde{G} & \tilde{F} & 0 & \cdots & 0 & 0 \\ 0 & \tilde{H} & \tilde{G} & \tilde{F} & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & \tilde{H} & \tilde{G} & \tilde{F} & 0 \\ 0 & 0 & \cdots & 0 & \tilde{H} & \tilde{G} & \tilde{F} \\ 0 & 0 & \cdots & 0 & 0 & \tilde{H} & \tilde{G} \end{bmatrix}_{(m-1)\times(m-1)} $$ $$ C= \begin{bmatrix} \tilde{H}u_0^{j}-Hu_0^{j+1} \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \\ \tilde{F}u_m^{j}-Fu_m^{j+1} \end{bmatrix}_{(m-1)\times 1} $$ $$ U^{j+1}= \begin{bmatrix} u_1^{j+1} \\ u_2^{j+1} \\ \vdots \\ u_{m-2}^{j+1} \\ u_{m-1}^{j+1} \end{bmatrix}_{(m-1)\times1} $$ $$ U^{j}= \begin{bmatrix} u_1^{j} \\ u_2^{j} \\ \vdots \\ u_{m-2}^{j} \\ u_{m-1}^{j} \end{bmatrix}_{(m-1)\times1} $$

Numerical Results

Example 1

Consider the diffusion equation obtained from Equation (1) by setting $\beta=\gamma=0$, with $\alpha=1$, $L=1$, and $t>0$. The exact solution for this example is given by

$$ u(x,t)=e^{-\pi^{2}t}\sin(\pi x) \tag{7} $$

The initial and boundary conditions are given by

$$ \left\{ \begin{array}{ll} u(x,0)=\sin(\pi x) & \text{Initial Condition} \\ u(0,t)=u(1,t)=0 & \text{Boundary Condition} \end{array} \right. \tag{8} $$

Fig: (a)Exact and Numerical Solution for Example-1 with different $h$ values $k= h^2$ (b) Behavior of Error for Example-1 with Different h values

Example 2

Consider the advection-diffusion equation obtained from Equation (1) by setting $\alpha = 0.1, \beta=1, \gamma=0$, with initial Condition given by

$$ \phi(x)=\exp(5x)\left[\cos\left(\frac{\pi}{2}x\right) +0.25\sin\left(\frac{\pi}{2}x\right)\right] \tag{9} $$

The the exact solution is given by

$$ u(x,t)= \exp\left(5\left(x-\frac{t}{2}\right)\right) \exp\left(-\frac{\pi^2}{40}t\right) \left[ \cos\left(\frac{\pi}{2}x\right) +0.25\sin\left(\frac{\pi}{2}x\right) \right] \tag{10} $$

The boundary conditions are obtained from the exact solution. To evaluate the accuracy of our proposed method,we solve this example for $L=1$ at the final time $T=2$

Fig: Numerical and Exact Solution for Example-2 with different $h$ values $k= h^2$

Example 3

Consider the advection-diffusion equation obtained from Equation (1) by setting $\alpha = \frac{\beta}{P_e}, \beta=-11, \gamma=0$, with different $P_e = 10,20,100,1000 $ alogn with initial Condition given by

$$ u(x,t)=a\exp(bt-\bar{c}x). \tag{9} $$

where

$$ \bar{c} = \frac{-\beta+\sqrt{\beta^{2}+4\alpha b}} {2\alpha}. \tag{12} $$

and the boundary conditions are given as follows:

$$ u(0,t)=a\exp(bt), \qquad u(1,t)=a\exp(bt-\bar{c}). \tag{13} $$

Fig: Numerical and Exact Solution for Example-3 with different $P_e$ values $k= h^2$

Example 4

Consider the advection-diffusion-reaction equation as given below

$$ \frac{\partial u}{\partial t}(x,t) + \frac{\partial u}{\partial x}(x,t) - 3u(x,t) = \frac{\partial^2 u}{\partial x^2}(x,t), \qquad (x,t)\in [0,2\pi]\times\left[0,\frac{1}{2}\right]. \tag{14} $$

with the initial condition

$$ u(x,0)=\sin(x). \tag{15} $$

Boundary condition

$$ u(0,t)=e^{2t}\sin(-t), \qquad u(2\pi,t)=e^{2t}\sin(2\pi-t). \tag{16} $$

Fig: Numerical and Exact Solution for Example-4 with different $h$ values $k= h^2$

Compact finite difference schemes for solving partial differential equations

Supervisor : Dr. Harshita Madduri ( NIT Kurukshetra)