The default Crank-Nicholson method (Implicit Scheme)

This post is part of a series of articles on the Finite Difference Method, with other posts in the series focusing on:

• Approximation of derivatives using the finite difference method
• Solving explicit diffusion equations
• Default Crank-Nicholson method
• Solving a tridiagonal matrix system using Thomas's algorithm

In the previous article about the Finite Difference Method, it was shown that the explicit method for numerically solving the heat equation requires a very small time step size, which prompts the consideration of a new method that can use larger time steps but at the cost of increased computational work per step. This method is called the Crank-Nicholson Scheme. 

Explicit method for the heat equation

The explicit method uses the forward difference for the time derivative and the centered second difference for the second spatial derivative:

uin+1​−uin/Δt​​=α - ui+1n​−2uin​+ui−1n​​2/(Δx)2

Implicit Scheme

The Crank-Nicholson Scheme improves this method by using a weighted average of the second-order derivative with respect to position at time n and n+1. The question that arises is: the value of uuu at time n+1 is still unknown. How can it be determined? This question will be answered below.

Another question is how to divide the weights of the second derivatives at both time intervals. The Crank-Nicholson method Use a 50-50 average, but other options are also available. However, the stability of the method is related to the weighting parameter θ in the general case:

θ=0⇒Explicit

θ=1⇒Fully Implicit

θ=0.5⇒Crank-Nicholson (midpoint)

The Crank-Nicholson scheme for the 1D heat equation

The Crank-Nicholson equation for the one-dimensional heat equation is as follows:

Δtuin+1​−uin​​=α21​((Δx)2ui+1n+1​−2uin+1​+ui−1n+1​​+(Δx)2ui+1n​−2uin​+ui−1n​​)

If it is defined that

λ=(Δx)²αΔt

This equation can be rearranged to separate the known term (current time) and the unknown term (next time) as follows:

−2λ​ui−1n+1​+(1+λ)uin+1​−2λ​ui+1n+1​=2λ​ui−1n​+(1−λ)uin​+2λ​ui+1n​

Since the left side has three unknown terms (next time), it means that at each time step, multiple linear equations must be solved simultaneously. For each time step nnn, there will be N−2 equations (if there are a total of N grid points) corresponding to the interior points of the grid area. By using Dirichlet boundary conditions (constant values), the endpoint points do not need to be solved in this system.

Linear System

A linear system can be represented in matrix form as:

Aun+1=Bun

Where:

• A and B are tridiagonal matrices.
• u^n and u^n+1 are vectors of temperature values at points inside the rod at times n and n+1.
• iii Starting from 1 to N−2

The structure of matrices A and B will be:

A=​1+λ−2λ​000​−2λ​1+λ⋱⋯⋯ 0−2λ​⋱−2λ​0​⋯⋯⋱1+λ−2λ​​000−2λ​1+λ​​

B will have the same structure, but replace 1+λ with 1−λ1 and −λ/2 with λ/2.

Solving a system of linear equations

The first attempt to solve this system might be to invert matrix A to find:

u^n+1=A^−1 Bu^n

But this approach is not suitable because it does not take into account the sparsity and band structure of matrix A. In the next article, we will see an efficient method for solving this system by utilizing its structure, so that the computation per time step uses the least resources.

Reference: Crank-Nicholson Implicit Scheme

From https://www.quantstart.com/articles/Crank-Nicholson-Implicit-Scheme/