Back to Projects

Analysis of methods to solve of Heat Equations

Supervisor : Dr. Harshita Madduri ( NIT Kurukshetra)

The following is a brief summary of the methodology and results of my Masters summer internship at NIT Kurukshetra.

Project Abstract

The numerical solution of the Heat Equation plays a crucial role in various scientific and engineering applications. In this study, we implemented Explicit, Implicit, and Crank-Nicholson schemes to solve the Heat Equation and performed a comparative analysis of their computational costs. To efficiently solve the resulting Tri-Diagonal System of Linear Equations (SOLE), we employed the Thomas Algorithm and Sparse Matrix methods, ensuring optimal computational efficiency. The implementation was carried out using Python, with a focus on accuracy, stability, and performance. The comparative results highlight the trade-offs between different schemes in terms of convergence, stability, and computational complexity, providing valuable insights for selecting the most appropriate method for practical applications.

Key Features

Numerical Solution of the One-Dimensional Heat Equation

In this section, we present the numerical solution of the one-dimensional heat equation. The heat equation is a fundamental partial differential equation that describes the distribution of temperature in a medium over time. We consider a model with specified initial and boundary conditions.

Mathematical Model

The governing heat equation is given by:

\[ \frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2} \]

where:

Boundary and Initial Conditions

The problem is solved subject to the following conditions:

\[ \begin{cases} u(0,t)=0, & t > 0, \\ u(L,t)=0, & t > 0, \end{cases} \]

These are homogeneous Dirichlet boundary conditions, indicating that the temperature at both ends of the rod remains fixed at zero.

\[ u(x,0)=f(x), \qquad 0 \le x \le L \]

This represents the initial temperature distribution along the domain.

Domain of the Problem

\[ 0 < x < L, \qquad t > 0 \]

The objective is to determine the temperature distribution \(u(x,t)\) over the spatial domain and time interval using a numerical method such as the Finite Difference Method (FDM).

Finite Difference Methods for the One-Dimensional Heat Equation

Consider the one-dimensional heat equation

\[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \]

where \(u(x,t)\) denotes the temperature distribution and \(\alpha\) is the thermal diffusivity coefficient.

\[ x_i=i\Delta x,\qquad t_j=j\Delta t \]

Let

\[ u(x_i,t_j)=u_i^j \]

Finite Difference Approximations

Forward Difference in Time

\[ \frac{\partial u}{\partial t} \approx \frac{u_i^{j+1}-u_i^j}{\Delta t} \]

Backward Difference in Time

\[ \frac{\partial u}{\partial t} \approx \frac{u_i^j-u_i^{j-1}}{\Delta t} \]

Central Difference in Space

\[ \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1}^{j}-2u_i^{j}+u_{i-1}^{j}} {(\Delta x)^2} \]
\[ \lambda= \frac{\alpha\Delta t} {(\Delta x)^2} \]

1. Explicit Finite Difference Scheme (FTCS)

Derivation

\[ \frac{u_i^{j+1}-u_i^j}{\Delta t} = \alpha \frac{u_{i+1}^{j}-2u_i^{j}+u_{i-1}^{j}} {(\Delta x)^2} \]
\[ u_i^{j+1}-u_i^j = \lambda \left( u_{i+1}^{j} -2u_i^{j} +u_{i-1}^{j} \right) \]
\[ u_i^{j+1} = \lambda u_{i-1}^{j} +(1-2\lambda)u_i^{j} +\lambda u_{i+1}^{j} \]

Stability Condition

\[ \lambda \le \frac{1}{2} \]
\[ \frac{\alpha\Delta t} {(\Delta x)^2} \le \frac12 \]

The Explicit FTCS method is conditionally stable and requires small time steps for convergence.

2. Bender–Schmidt Method

The Bender–Schmidt method is a special case of the Explicit scheme obtained by choosing

\[ \lambda=\frac12 \]

Substituting into the explicit formula:

\[ u_i^{j+1} = \frac12u_{i-1}^{j} + \frac12u_{i+1}^{j} \]
\[ u_i^{j+1} = \frac12 \left( u_{i-1}^{j} + u_{i+1}^{j} \right) \]

Characteristics

Stability Condition

\[ \lambda=\frac12 \]

The method remains stable only when this condition is satisfied.

3. Implicit Finite Difference Scheme (BTCS / Laasonen Method)

Derivation

\[ \frac{u_i^{j+1}-u_i^j} {\Delta t} = \alpha \frac{ u_{i+1}^{j+1} -2u_i^{j+1} +u_{i-1}^{j+1} } {(\Delta x)^2} \]
\[ u_i^{j+1}-u_i^j = \lambda \left( u_{i+1}^{j+1} -2u_i^{j+1} +u_{i-1}^{j+1} \right) \]
\[ -\lambda u_{i-1}^{j+1} +(1+2\lambda)u_i^{j+1} -\lambda u_{i+1}^{j+1} = u_i^j \]

Matrix Form

\[ AU^{j+1}=U^j \]

Stability

\[ \text{BTCS is Unconditionally Stable} \]

There is no restriction on the ratio \( \alpha\Delta t/(\Delta x)^2 \).

4. Crank–Nicolson Method

The Crank–Nicolson method is obtained by averaging the Explicit and Implicit finite difference schemes.

Derivation

\[ \frac{u_i^{j+1}-u_i^j} {\Delta t} = \frac{\alpha}{2} \left[ \frac{ u_{i+1}^{j} -2u_i^{j} +u_{i-1}^{j} } {(\Delta x)^2} + \frac{ u_{i+1}^{j+1} -2u_i^{j+1} +u_{i-1}^{j+1} } {(\Delta x)^2} \right] \]
\[ -\frac{\lambda}{2}u_{i-1}^{j+1} +(1+\lambda)u_i^{j+1} -\frac{\lambda}{2}u_{i+1}^{j+1} = \frac{\lambda}{2}u_{i-1}^{j} +(1-\lambda)u_i^{j} +\frac{\lambda}{2}u_{i+1}^{j} \]

Matrix Form

\[ AU^{j+1}=BU^j \]

Stability

\[ \text{Crank–Nicolson is Unconditionally Stable} \]

The method remains stable for all values of \( \lambda \).

Comparison of Numerical Schemes

Method Type Stability Condition Accuracy
Explicit (FTCS) Explicit Conditional \(\lambda \le 0.5\) O(Δt)+O(Δx²)
Bender–Schmidt Explicit Conditional \(\lambda=0.5\) O(Δt)+O(Δx²)
BTCS Implicit Unconditional No Restriction O(Δt)+O(Δx²)
Crank–Nicolson Semi-Implicit Unconditional No Restriction O(Δt²)+O(Δx²)

Sparse Matrix Solver

The coefficient matrices arising from the BTCS and Crank–Nicolson schemes are highly sparse because only the main diagonal and the two adjacent diagonals contain non-zero entries.

Instead of storing all matrix elements, sparse matrix techniques store only the non-zero values and their indices. This significantly reduces memory consumption and computational overhead.

In Python, sparse matrices can be efficiently represented using libraries such as SciPy:

\[ Ax=b \]

where the matrix \(A\) is stored in sparse format and solved using optimized sparse linear algebra routines.

Performance Analysis of Linear System Solvers

The Implicit (BTCS) and Crank–Nicolson finite difference schemes generate a system of linear equations at each time step. The resulting coefficient matrix is tridiagonal, making specialized numerical solvers highly effective. In this project, two approaches are considered:

  • Thomas Algorithm (Tri-Diagonal Matrix Algorithm)
  • Sparse Matrix Solver

1. Thomas Algorithm (Tri-Diagonal Matrix Algorithm)

Introduction

The Thomas Algorithm is a simplified form of Gaussian Elimination specifically designed for tridiagonal systems of linear equations. Since only three diagonals contain non-zero entries, the algorithm can solve the system efficiently using forward elimination and backward substitution.

General Tridiagonal System

The system obtained from the finite difference discretization can be written as:

\[ A\mathbf{u}=\mathbf{b} \]
\[ \begin{bmatrix} b_1 & c_1 & 0 & 0 & \cdots \\ a_2 & b_2 & c_2 & 0 & \cdots \\ 0 & a_3 & b_3 & c_3 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \\ 0 & 0 & \cdots & a_n & b_n \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_n \end{bmatrix} = \begin{bmatrix} d_1 \\ d_2 \\ d_3 \\ \vdots \\ d_n \end{bmatrix} \]

Equivalently,

\[ a_i u_{i-1} + b_i u_i + c_i u_{i+1} = d_i \]

Forward Elimination

\[ c'_1=\frac{c_1}{b_1} \]
\[ d'_1=\frac{d_1}{b_1} \]

For \(i=2,3,\ldots,n\):

\[ c'_i= \frac{c_i} {b_i-a_i c'_{i-1}} \]
\[ d'_i= \frac{d_i-a_i d'_{i-1}} {b_i-a_i c'_{i-1}} \]

This process removes the lower diagonal elements and transforms the matrix into an upper triangular form.

Backward Substitution

\[ u_n=d'_n \]
\[ u_i= d'_i-c'_i u_{i+1} \qquad (i=n-1,n-2,\ldots,1) \]

Computational Complexity

\[ \text{Time Complexity}=O(n) \]
\[ \text{Memory Complexity}=O(n) \]

Therefore, the Thomas Algorithm is significantly faster than general dense matrix solvers and is well suited for finite difference discretizations of the heat equation.

Application to BTCS and Crank–Nicolson Methods

For the BTCS scheme:

\[ -\lambda u_{i-1}^{j+1} +(1+2\lambda)u_i^{j+1} -\lambda u_{i+1}^{j+1} = u_i^j \]

Thus,

\[ a_i=-\lambda, \qquad b_i=1+2\lambda, \qquad c_i=-\lambda \]

which directly forms a tridiagonal system that can be solved using the Thomas Algorithm.

2. Sparse Matrix Solver

Introduction

In large-scale numerical simulations, the coefficient matrix contains a very small number of non-zero entries compared to the total number of elements. Such matrices are known as sparse matrices.

Instead of storing every element, sparse matrix methods store only the non-zero values and their locations, resulting in substantial savings in memory and computational cost.

Example of a Sparse Matrix

\[ A= \begin{bmatrix} 2 & -1 & 0 & 0 \\ -1 & 2 & -1 & 0 \\ 0 & -1 & 2 & -1 \\ 0 & 0 & -1 & 2 \end{bmatrix} \]

Only the main diagonal and the two adjacent diagonals contain non-zero values.

Sparse Representation

Rather than storing all matrix entries, only the following diagonals are stored:

\[ \text{Main Diagonal} = [2,2,2,2] \]
\[ \text{Upper Diagonal} = [-1,-1,-1] \]
\[ \text{Lower Diagonal} = [-1,-1,-1] \]

This significantly reduces storage requirements for large systems.

Mathematical Formulation

\[ A\mathbf{u} = \mathbf{b} \]

The matrix \(A\) is stored in sparse format while preserving the original mathematical problem.

Computational Complexity

\[ \text{Dense Solver} = O(n^3) \]
\[ \text{Dense Memory} = O(n^2) \]
\[ \text{Sparse Solver} \approx O(n) \text{ to } O(n\log n) \]
\[ \text{Sparse Memory} = O(n) \]

Sparse matrix methods become increasingly advantageous as the number of spatial grid points increases.

Relation to Python Scientific Computing

In Python, sparse matrices are commonly stored using CSR (Compressed Sparse Row) or CSC (Compressed Sparse Column) formats. Libraries such as SciPy provide optimized sparse linear algebra routines that solve

\[ A\mathbf{u} = \mathbf{b} \]

by operating only on the non-zero entries of the matrix. Internally, these methods employ sparse LU factorization and highly optimized numerical routines.

3. Performance Comparison

Both the Thomas Algorithm and Sparse Matrix Solver solve the same system of equations generated by BTCS and Crank–Nicolson methods. However, their implementations differ significantly.

Property Thomas Algorithm Sparse Matrix Solver
Matrix Type Tridiagonal Only General Sparse Matrix
Time Complexity O(n) O(n) – O(n log n)
Memory Usage O(n) O(n)
Flexibility Low High
Large Scale Problems Good Excellent
Parallelization Limited Supported

Observed Results

Numerical experiments conducted in this project showed that the Sparse Matrix Solver achieved approximately 10× faster execution time compared to the conventional Thomas Algorithm implementation for large grid sizes.

\[ T_{\text{sparse}} \ll T_{\text{Thomas}} \]

The improvement is mainly due to optimized sparse linear algebra routines, efficient memory access, reduced storage requirements, and highly optimized low-level implementations.

Conclusion

The Thomas Algorithm is an efficient direct solver for tridiagonal systems arising from finite difference discretizations. However, sparse matrix techniques provide greater flexibility and superior performance for large-scale simulations. In this project, the sparse matrix approach demonstrated significantly reduced execution times while maintaining the same numerical accuracy, making it the preferred choice for solving large heat equation problems.

Numerical Example

Consider the one-dimensional heat equation

\[ \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}, \qquad 0 < x < 1,\; t>0 \]

with homogeneous Dirichlet boundary conditions

\[ u(0,t)=0, \qquad u(1,t)=0, \qquad t>0 \]

and initial condition

\[ u(x,0) = 2\sin\left(\frac{\pi x}{2}\right), \qquad 0 \le x \le 1. \]
Mesh Points
N = M		 Error in Solution
Calculating by Thomas Algorithm		 Computational Cost
(using TA)		 Error in Solution
Calculating by Sparse Matrix		 Computational Cost
(Using SM)
10 0.5811 0.0019919872283933547 0.4344 0.0009987354278564453
20 0.3915 0.02094101905822754 0.2651 0.0009996891021728516
30 0.2390 0.06778788566589355 0.1619 0.0009913444519042969
40 0.14594 0.16256427764892578 0.0988 0.00099135485609412969
Example-4
Github Repository
Back to Projects