[Python] Fluid simulation: Implement the diffusion equation

Introduction

In addition to studying computational fluid dynamics (CFD), which is a study related to simulation of fluids such as air and water, I would like to summarize the knowledge necessary for constructing computational fluid dynamics (CFD) for water (in multiple articles).

I feel that computational fluid dynamics is a rather difficult discipline, so I would like to write it in a way that is easy for beginners to understand. It seems that there are many mistakes, so please contact us if you find it. Also, I would appreciate it if you could comment on what is difficult to understand. We will update it from time to time.

Target audience

People who can use Python
Those who are interested in numerical calculation
Those who are interested in fluid mechanics
Those who understand basic university physics and mathematics (differential equations?)

Those who are interested in matrix calculation

series

Chapter 1: [Python] Fluid Simulation: Implementing Advection Equation
Chapter 2: This article
Supplementary position in Chapter 2: [Python] Article that enables high-speed sparse matrix calculation
Chapter 3: [Python] Fluid Simulation: From Linear to Nonlinear
Chapter 4: [Python] Fluid Simulation: Incompressible Navier-Stokes Equation

Rough content of this article

As a preliminary step to deal with the basic fluid equations required for water simulation, we will also briefly summarize and implement the diffusion equations. ** I've made it understandable without reading the previous article (probably) **

Table of contents in the article
[Diffusion equation](#Diffusion equation (formula for uniform field))
[Diffusion number-based stability judgment](#Diffusion number-based stability judgment)
[Calculation conditions](# Calculation conditions)
Implementation

[Center difference](# 1-Center difference)
[Crank Nicholson's implicit method](# 2-Crank Nicholson's implicit method)
[Direct method](# 2-1-Direct method)
[Iterative method](# 2-2-Iterative method)
[Standing Iterative Method](# 2-2-1-Standing Iterative Method)
[Unsteady Iterative Method](# 2-2-2-Unsteady Iterative Method krylov Subspace Method)
[Implementation of stationary iterative method](# 2-3-Implementation of stationary iterative method)
[Jacobi method](# 2-3-1-jacobi method)
[Simple example](# 2-3-1-1-Simple example)
[Caution](# 2-3-1-2-jacobi method precautions)
[Implementation](# 2-3-1-3-jacobi method implementation)
[Gauss-Seidel method](# 2-3-2-gauss-seidel method Gauss-Seidel method)
[Caution](# 2-3-2-1-gauss-seidel method precautions)
[Implementation](# 2-3-2-2-2-gauss-seidel method implementation)
[SOR method](# 2-3-3-sor method)
[Implementation](# 2-3-3-1-sor method implementation)

Summary
References

Diffusion equation (equation that expresses the uniformization of the field)

What is the one-dimensional diffusion equation? $ \frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2} $

It is expressed by the formula. In terms of physical meaning, it means that physical quantities are uniformly uniform. It means heat conduction and diffusion of substances as well as heat. For example, when you put milk in coffee, it spreads slowly without stirring, or when you put the end of a metal stick in a boiling pot as shown in the figure below, the whole stick gradually becomes hot water. The phenomenon of the same temperature corresponds to diffusion.

Example of heat diffusion
The upper one-dimensional diffusion equation is an equation that predicts the temperature distribution of the metal rod in the situation shown in the figure below.
$ \ Kappa $ in the formula changes the speed at which the temperature propagates. The larger it is, the shorter it takes for the entire bar to reach 100 ° C.

Methods for discretizing this one-dimensional diffusion equation include the explicit difference method, the central finite difference method, and the implicit method, Crank Nicholson's implicit method. As a reminder, the explicit method predicts the future value one hour later based on the known value of the current time, and the implicit method predicts the future value one hour later as well. We will implement these two. If you write down each with the difference formula, it will be as follows.

Center difference It can be derived by performing the first-order accuracy upwind difference explained in the previous article twice. It is a kind of explicit method that predicts the value of the next time n + 1 only with the data of the current time n. $ \frac{T_j^{n+1} - T_j^n}{\Delta t} = \kappa \frac{T_{j+1}^n - 2 T_j^n + T_{j-1}^n }{\Delta x^2} $
Crank Nicholson's implicit method

A method of averaging the partial differential with respect to the spatial derivative between the known value (time n) and the unknown value one step later (time n + 1). As the name implies, it is an implicit method that also uses the value of the next time n + 1 to predict the value of the next time n + 1. Details will be described later. $ \frac{T_j^{n+1} - T_j^n}{\Delta t} = \frac{\kappa}{2} \left(\frac{T_{j+1}^n - 2 T_j^n + T_{j-1}^n }{\Delta x^2} + \frac{T_{j+1}^{n+1} - 2 T_j^{n+1} + T_{j-1}^{n+1} }{\Delta x^2}\right) $

Determining stability by diffusion number

Of the discretization methods shown above, the ** explicit method using the central difference ** $ d = \kappa \frac{\Delta t}{\Delta x^2} $ It is necessary to determine the stability of the numerical calculation by using what is called the ** diffusion number d **. On the other hand, the implicit method is unconditionally stable.

In particular,

d \leq \frac{1}{2} \quad or \quad \Delta t \leq \frac{(\Delta x)^2}{2 \kappa}

Must be within the range of. The important thing in this ** diffusion number condition is that if you want to reduce the step size $ \ Delta x $, you must reduce the time step size $ \ Delta t $ by its square **. For example, if you multiply $ \ Delta x $ by 1/2, you must reduce $ \ Delta t $ to 1/4. As a result, solving the diffusion equation by the explicit method tends to increase the computational load, so the implicit method is generally used in many cases.

In addition, this condition is obtained from the stability analysis of Von Neumann, but I will not explain it here, but only an intuitive image.

If you rewrite the formula of the central difference,

T_j^{n+1} = d T_{j+1}^n + (1 - 2d) T_j^n + d T_{j-1}^n

It will be. Diffusion means that it is uniformed in a disorderly manner, so when the coefficient $ (1-2d) $ of $ T_j ^ n $ becomes negative, it passes more than the physical quantity of the jth lattice to the next. Therefore, it should be $ d \ leq0.5 $.

I think it's hard to understand, so I'll explain it with the example of money below (** It's just an example to have an image ). Suppose you have a group of N people, and at a certain time n, the j-1st person has 30 yen, the jth person has 100 yen, and the j + 1st person has 50 yen. Diffusion is an operation that allows these N people to have the same amount of money ( uniformization **), so keep exchanging a certain percentage of money with neighbors until the same amount of money is reached. I will. If you look closely at the formula for the central difference, this constant percentage is the diffusion number d. Therefore, if you set d = 0.1 and calculate, the jth person will have 88 yen at the n + 1 time one hour later. Given this, it's easy to see that the diffusion number d must not be greater than 1/2.

** * This is an example to grasp the atmosphere of diffusion **

Calculation condition

Problem
Consider a substance (such as metal) with a temperature gradient of 100K to 200K, and calculate the temperature change when a heat source of 150K is applied to both ends. As you can easily predict, the whole substance will be constant at 150K after enough time. Let's predict the temperature distribution in the middle of the change using the one-dimensional thermal diffusion equation. (Supplement: K is a unit of temperature called Kelvin, which corresponds to 273.15K = 0 ° C)

The lattice width $ \ Delta x $, thermal conductivity $ \ kappa $, and time step width $ \ Delta t $ are all calculated as constants. This time, $ \ Delta x = 1 $, $ \ Delta t = 0.2 $, $ \ kappa = 0.5 $, and the diffusion number $ d = 0.1 $ is calculated under stable conditions.

import numpy as np
import matplotlib.pyplot as plt
Num_stencil_x = 101
x_array = np.float64(np.arange(Num_stencil_x))
temperature_array = x_array + 100
temperature_lower_boundary = 150
temperature_upper_boundary = 150
Time_step = 100
Delta_x = max(x_array) / (Num_stencil_x-1)
Delta_t = 0.2
kappa = 0.5
d = kappa * Delta_t / Delta_x**2
exact_temperature_array = (temperature_upper_boundary - temperature_lower_boundary) / (x_array[-1] - x_array[0]) * x_array + temperature_lower_boundary
plt.plot(x_array, temperature_array, label="Initial condition")
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("Temperature [K]")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

Implementation

1. Center difference

\frac{T_j^{n+1} - T_j^n}{\Delta t} = \kappa \frac{T_{j+1}^n - 2 T_j^n + T_{j-1}^n }{\Delta x^2}

The figure below shows the result when the diffusion number d is 0.1. I tried to output the graph by dividing it into several time steps (time). After 100 hours, both ends have reached almost 150K, and the temperature distribution is dragged by it like a sine curve. After 5000 hours, the temperature of the entire bar is asymptotic to 150K, indicating that the implementation of the center difference is successful. Since the implementation is easy to follow expressions, I will write it without using numpy functions as much as possible.

By the way, if you set $ \ kappa = 5 $ and calculate with the diffusion number $ d = 1 $, you can see that the calculation diverges and the solution cannot be found as shown in the figure below.


temperature_explicit = temperature_array.copy()
for n in range(Time_step):
    temperature_old = temperature_explicit.copy()
    temperature_explicit[0] += kappa * Delta_t / Delta_x**2 * \
        (temperature_explicit[1] - 2*temperature_old[0] + temperature_lower_boundary)
    temperature_explicit[-1] += kappa * Delta_t / Delta_x**2 * \
        (temperature_upper_boundary - 2*temperature_old[-1] + temperature_old[-2])
    for j in range(1, Num_stencil_x-1):
        temperature_explicit[j] += kappa * Delta_t / Delta_x**2 * \
            (temperature_old[j+1] - 2*temperature_old[j] + temperature_old[j-1])
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.plot(x_array, temperature_explicit, label="Explicit(100 steps)")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

2. Crank Nicholson's implicit method

\frac{T_j^{n+1} - T_j^n}{\Delta t} = \frac{\kappa}{2} \left(\frac{T_{j+1}^n - 2 T_j^n + T_{j-1}^n }{\Delta x^2} + \frac{T_{j+1}^{n+1} - 2 T_j^{n+1} + T_{j-1}^{n+1} }{\Delta x^2}\right)

If you transform this and bring the value of time n + 1 to the left side $ -T_{j+1}^{n+1} +2 \left(\frac{1}{d} + 1 \right) T_j^{n+1} - T_{j-1}^{n+1} = T_{j+1}^{n} +2 \left(\frac{1}{d} - 1 \right) T_j^{n} + T_{j-1}^{n} \\ d = \kappa \frac{\Delta t}{\Delta x^2} $

Assuming that the grid points exist in the range of 1 to M, the boundary value is represented by $ T_0, T_ {M + 1} $, and it is converted to a matrix display.


\left(
\begin{array}{cccc}
      2 \left(\frac{1}{d} + 1 \right) & -1 &  0 & \ldots  & \ldots & 0 \\
      -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots & 0 \\
      0  &-1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots  \\
      \vdots & \ddots & \ddots & \ddots & \ddots & \ddots\\
      0  & 0 & \ddots & -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 \\
      0  & 0 & \ldots & 0 & -1 & 2 \left(\frac{1}{d} + 1 \right)
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      T_1^{n+1}  \\
      T_2^{n+1}  \\
      T_3^{n+1}    \\
      \vdots \\
      T_{M-1}^{n+1} \\
      T_M^{n+1} 
\end{array}
\right)


= \left( \begin{array}{c}
    T_2^{n} + 2 \left(\frac{1}{d} - 1 \right) T_1^{n} + \left(T_0^n + T_0^{n+1}\right) \\
      T_3^{n} + 2 \left(\frac{1}{d} - 1 \right) T_2^{n} + T_1^n  \\
      T_4^{n} + 2 \left(\frac{1}{d} - 1 \right) T_3^{n} + T_2^n  \\
      \vdots \\
      T_M^{n} + 2 \left(\frac{1}{d} - 1 \right) T_{M-1}^{n} + T_{M-2}^n  \\
      \left(T_{M+1}^n + T_{M+1}^{n+1}\right) + 2 \left(\frac{1}{d} - 1 \right) T_{M}^{n} + T_{M-1}^n
    \end{array} \right)

Consider solving this simultaneous one-dimensional equation. When the number of grid points M is small, the unknown is small, so a direct method that directly finds a solution, such as a method called Gauss's direct method, is used. However, the larger M is, the more difficult it is to find a solution in terms of computational complexity. Therefore, in general, in order to solve a large-scale simultaneous one-dimensional equation, a method called ** iterative method ** is used to converge the approximate solution to an exact solution.

2-1. Direct method

The types are as follows. I will explain when it is needed.

Gaussian elimination
The idea is the same as the sweeping method learned in the linear algebra class at the university.
LU decomposition
Cholesky decomposition
Cholesky decomposition
Modified Cholosky decomposition
Incomplete Cholosky decomposition (IC method)

2-2. Iterative method

Ax=b

Consider solving the one-dimensional simultaneous equations. this $ r=b-Ax' $ The method of finding an approximate solution of the exact solution $ x ^ * = A ^ {-1} b $ is repeated by iteratively changing the estimated value $ x'$ until r becomes sufficiently small. It is a law. The iterative method can be divided into the stationary method and the non-stationary method (Krylov subspace method).

2-2-1. Steady Iterative Method

The stationary method (steady iterative method) is a method that does not change except for the solution vector x during iterative calculation. Since it is generally slow, either use the non-stationary method described later when performing numerical calculations, or use it together as a pre-processing of the non-stationary method (processing performed to find an approximate approximate solution). think.

The following three methods can be mentioned as typical methods. In the following, the estimated value of the solution in m iterations is expressed as $ x ^ {(m)} $, and the estimated value $ x ^ {(m)} $ in the mth iteration is known as the m + 1th iteration. Here is an example of the calculation to find the estimated value of $ x ^ {(m + 1)} $.

Jacobi method

A method of calculating the estimated solution $ x_i ^ {(m + 1)} $ in the i-th row of the m + 1 iteration using only the known estimated solution $ x_i ^ {(m)} . $ x_i^{(m+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j=1,j\not=i}^{n} a_{ij} x_j^{(m)} \right) $$ 2. Gauss-Seidel method

M + 1 already calculated as the known estimated solution $ x_i ^ {(m)} $ when finding the estimated solution $ x_i ^ {(m + 1)} $ on the i-th row of the m + 1 iteration. A method of calculation using the estimated solution $ x_ {0} ^ {(m + 1)}, \ cdots, x_ {i-1} ^ {(m + 1)} $ of the iteration. Basically, the convergence calculation finishes faster than the Jacobi method. $ x_i^{(m+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(m+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(m)} \right) $ 3. SOR (Successive Over-Relaxation) method / Sequential acceleration relaxation method

Gauss-Seidel method plus relaxation factor $ \ omega $. When the relaxation coefficient is $ \ omega = 1 , the Gauss-Seidel method is used. Also, if it is 1 or less, convergence is slower than Gauss-Seidel, but problems that cannot be solved by the Gauss-Seidel method can be solved. Basically, set it to a value of 1 or more, but if it is set too large, it will diverge, so it is important to select an appropriate value. It seems that 1.9 is often selected. $ x_i^{(m+1)} = x_i^{(m)} + \omega \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(m+1)} - \sum_{j=i}^{n} a_{ij} x_j^{(m)} \right) $$

Multigrid method There are the geometric multi-grid method and the algebraic multi-grid (AMG) method. The latter AMG method is often used as a pretreatment for the transient method. I think that many people will know it for the first time, so I will write the details in another article.

Details of the former three items will be described later.

2-2-2. Transient Iterative Method (Krylov Subspace Method)

Is it like a good boy with direct and iterative methods? It is classified as an iterative method.

The most famous non-stationary method is the Conjugate Gradient method (CG method). Detailed explanation is given in this article. There is a figure and it is insanely easy to understand. In addition, this article explains the difference from the steepest descent method that is often used in machine learning. The image of the CG method is also written and it is easy to understand.

Symmetric matrix
Conjugate Gradient method (CG method)
Preconditioned Conjugate Gradient method (PCG method)
Incomplete Cholesky Conjugate Gradient method (ICCG method)
Incomplete LU Conjugate Gradient method (ILUCG)
Other
Asymmetric matrix
Bi-Conjugate Gradient (BiCG method)
Conjugate Residual Squared (CGS method)
BiCG Stabilization (BiCG STabilization)
Generalized Conjugate Residual (GCR method)
Generalized Minimal Residual (GMRES method)

These techniques are basically provided as a library. Also, I think there are a lot of calculation codes on Github. If you want to study in more detail, please refer to that as well. In another article, I would like to explain and implement as many methods as possible for the transient method.

For the time being, in this article, I would like to use the Jacobi method, Gauss-Seidel method, and SOR method, which are stationary iterative methods, among the algorithms explained so far.

2-3. Implementation of stationary iterative method

The implementation will be carried out with reference to the Material given by Professor Nakajima of the University of Tokyo. If you find the explanation difficult to understand, we would appreciate it if you could refer to it.

In the example below

A = \left(
    \begin{array}{ccc}
      a_{1,1} & a_{1,2} & \ldots &a_{1,n}  \\
      a_{2,1} & a_{2,2} & \ldots &a_{2,n} \\
      \vdots & \vdots &  \ddots &\vdots \\
       a_{n,1} & \ldots &  \ldots &a_{n,n}
    \end{array}
  \right) \\
D = \left(
    \begin{array}{ccc}
      a_{1,1} & 0 & \ldots &0  \\
      0 & a_{2,2} & \ldots &0 \\
      \vdots & \vdots &  \ddots &\vdots \\
       0 & \ldots &  \ldots &a_{n,n}
    \end{array}
  \right) \\
L = \left(
    \begin{array}{ccc}
      0 & 0 & \ldots &0  \\
      a_{2,1} & 0 & \ldots &0 \\
      \vdots & \vdots &  \ddots &\vdots \\
       a_{n,1} & \ldots &  \ldots &0
    \end{array}
  \right) \\
U = \left(
    \begin{array}{ccc}
      0 & a_{1,2} & \ldots &a_{1,n}  \\
      0 & 0 & \ldots &a_{2,n} \\
      \vdots & \vdots &  \ddots &\vdots \\
       0 & \ldots &  \ldots &0
    \end{array}
  \right)

Define the matrix as follows. D is only the diagonal component of A, L is the lower side of the diagonal component, and U is the upper side of the diagonal component.

2-3-1. Jacobi method

When $ A = D + L + U $ using the matrix D with the diagonal components of the matrix A, the upper matrix U, and the lower matrix L,

　Ax = b\\
 (D+L+U) x = b\\
 Dx = b - (L+U)x \\
 x = D^{-1}(b-(L+U)x)

$ Ax = b $ can be transformed like this, and the method of iteratively finding an approximate solution with the formula at the bottom is called the Jacobi method.

x_i^{(m+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j=1,j\not=i}^{n} a_{ij} x_j^{(m)} \right)

2-3-1-1. Simple example

Let A be a $ 3 \ times3 $ matrix


\left(
    \begin{array}{ccc}
      a_{1,1} & a_{1,2} &  a_{1,3}  \\
      a_{2,1} & a_{2,2} &  a_{2,3} \\
      a_{3,1} & a_{3,2} &  a_{3,3}
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      x_{1} \\
      x_{2} \\
      x_{3} 
    \end{array}
  \right)
  =
  \left(
    \begin{array}{c}
      b_{1} \\
      b_{2} \\
      b_{3} 
    \end{array}
  \right)

Consider solving the one-dimensional simultaneous equations.

Expressing the number of iterations as (m), the value of $ x ^ {(m + 1)} $ at (m + 1) after one iteration is

x_1^{(m+1)} = \left(b_1 - a_{1,2} x_2^{(m)} - a_{1,3} x_3^{(m)} \right) / a_{1,1} \\
x_2^{(m+1)} = \left(b_2 - a_{2,1} x_1^{(m)} - a_{2,3} x_3^{(m)} \right) / a_{2,2} \\
x_3^{(m+1)} = \left(b_3 - a_{3,1} x_1^{(m)} - a_{3,2} x_2^{(m)} \right) / a_{3,3}

It can be expressed in the form of. As you can see from the form of this equation, $ x ^ {(m + 1)} = D ^ {-1} (b-(L + R) x ^ {(m)}) $ holds. I will.

Then, by repeating the above calculation until the convergence condition shown below is satisfied, the exact solution $ x ^ * $ can be obtained.

\sum_{j=1}^{n} \left| \frac{x_j^{(m+1)}-x_j^{(m)}}{x_j^{(m+1)}} \right|= \epsilon

However, $ \ epsilon $ should be a sufficiently small positive number.

Now, let's practice with the Jacobi method with a simple problem.

As a simple example

\left(
    \begin{array}{ccc}
      3 & 2 &  -0.5  \\
      1 & 4 &  1  \\
      -1 & 0 &  4  
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      x_1   \\
      x_2  \\
      x_3
    \end{array}
  \right)
  =
  \left(
    \begin{array}{c}
      3   \\
      2  \\
      1
    \end{array}
  \right)

Think about solving. The answer is

\left(
    \begin{array}{c}
      1   \\
      0.125  \\
      0.5
    \end{array}
  \right)

Unlike before, we will implement it by making heavy use of numpy functions. Because it is overwhelmingly easy to see. The implementation is as follows. When you actually calculate

Solution [1.00000063 0.12499957 0.50000029]
Answer [1.    0.125 0.5  ]

Is displayed and it can be confirmed that the answer is being sought.

def jacobi(a_matrix, b_array, target_residual):
    """
Jacobi method
    Ax = b
    A = (D+L+U)Then
    Dx = b - (L+U)x
    x = D^{-1} (b - (L+U)x)
However, D is a diagonal matrix, L+Let U be the residual matrix.
    
    Parameters
    ----------
    a_matrix : numpy.float64
m × n matrix
    b_array : numpy.float64
m-row matrix
    target_residual : numpy.float64
Positive number. Target residuals.

    Returns
    -------
    x : numpy.float64
m-row matrix
    """
    x = b_array.copy()
    x_old = b_array.copy()
    
    diag_matrix= np.diag(a_matrix)  #Diagonal matrix
    l_u_matrix = a_matrix - np.diagflat(diag_matrix)  #Residual matrix
    count = 0
    while True:
        count += 1
        x = (b_array - np.dot(l_u_matrix, x_old))/diag_matrix
        residual = np.linalg.norm(x - x_old) / np.linalg.norm(x)
        x_old = x
        if residual <= target_residual:
            break
        elif count >= 10000:
            import sys
            print(residual)
            sys.exit()
    return x

A = np.array([[ 3.0,  2.0, -0.5],
              [ 1.0,  4.0,  1.0],
              [-1.0,  0.0,  4.0]])
b = np.array([3.0, 2.0, 1.0])

x = jacobi(A, b, 10**-6)
print("Solution", x)
print("Answer", np.dot(np.linalg.inv(A),b))

2-3-1-2. Precautions for Jacobi method

The Jacobi method (and Gauss-Seidel method) does not converge unless the matrix A satisfies the diagonal advantage (in the narrow sense). To briefly state what the diagonal dominance (in the narrow sense) is, we say that the diagonal component is greater than the sum of the absolute values of the other components in the same row. In the exercise above, try playing with the non-diagonal components larger.

\left| a_{ii} \right| > \sum_{j=1, j\not=i}^{n} \left| a_{ij} \right|

2-3-1-3. Jacobi method implementation

I will repost the formula to be solved.


\left(
\begin{array}{cccc}
      2 \left(\frac{1}{d} + 1 \right) & -1 &  0 & \ldots  & \ldots & 0 \\
      -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots & 0 \\
      0  &-1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots  \\
      \vdots & \ddots & \ddots & \ddots & \ddots & \ddots\\
      0  & 0 & \ddots & -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 \\
      0  & 0 & \ldots & 0 & -1 & 2 \left(\frac{1}{d} + 1 \right)
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      T_1^{n+1}  \\
      T_2^{n+1}  \\
      T_3^{n+1}    \\
      \vdots \\
      T_{M-1}^{n+1} \\
      T_M^{n+1} 
\end{array}
\right)


= \left( \begin{array}{c}
    T_2^{n} + 2 \left(\frac{1}{d} - 1 \right) T_1^{n} + \left(T_0^n + T_0^{n+1}\right) \\
      T_3^{n} + 2 \left(\frac{1}{d} - 1 \right) T_2^{n} + T_1^n  \\
      T_4^{n} + 2 \left(\frac{1}{d} - 1 \right) T_3^{n} + T_2^n  \\
      \vdots \\
      T_M^{n} + 2 \left(\frac{1}{d} - 1 \right) T_{M-1}^{n} + T_{M-2}^n  \\
      \left(T_{M+1}^n + T_{M+1}^{n+1}\right) + 2 \left(\frac{1}{d} - 1 \right) T_{M}^{n} + T_{M-1}^n
    \end{array} \right)

The result itself is about the same as the central difference. The diffusion number d is 0.1.

The figure below shows the result of calculation with the diffusion number d set to 1 as in the case of the center difference. ** You can see that the solution can be found firmly even if the diffusion number d is 0.5 or more **. This is the advantage of using the implicit method. Also, from a physical point of view, as the diffusion number d increases, the time step width and the space step width are constant, so the heat transfer speed becomes faster and converges to 150K faster than when the diffusion number d = 0.1. You can see that it is doing.

An implementation example is shown below. The jacobi function in the syntax was built in the previous section.


temperature_jacobi = temperature_array.copy()
for n in range(Time_step):
    a_matrix = np.identity(len(temperature_jacobi)) * 2 *(1/d+1) \
                - np.eye(len(temperature_jacobi), k=1) \
                - np.eye(len(temperature_jacobi), k=-1)
    temp_temperature_array = np.append(np.append(
                        temperature_lower_boundary, 
                        temperature_jacobi), temperature_upper_boundary)
    b_array = 2 * (1/d - 1) * temperature_jacobi + temp_temperature_array[2:] + temp_temperature_array[:-2]
    b_array[0] += temperature_lower_boundary
    b_array[-1] += temperature_upper_boundary
    temperature_jacobi = jacobi(a_matrix, b_array, 1e-8)
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.plot(x_array, temperature_jacobi, label="Implicit Jacobi(100 steps)")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

2-3-2. Gauss-Seidel method (Gauss-Seidel method)

When $ A = D + L + U $ using the matrix D having diagonal components of the matrix A, the upper triangular matrix U, and the lower triangular matrix L,

　Ax = b\\
 (D+L+U) x = b\\
 (D+L)x^{(m+1)} = b - Ux^{(m)} \\
 x^{(m+1)} = D^{-1}(b- L x^{(m+1)} - U x^{(m)})\\
x^{(m+1)} = (D+L)^{-1}(b - U x^{(m)})

The Gauss-Seidel method is a method that can transform $ Ax = b $ and find an approximate solution iteratively with the formula at the bottom.

x_i^{(m+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(m+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(m)} \right)

2-3-2-1. Gauss-Seidel method precautions

As mentioned in the Jacobi method, the Gauss-Seidel method does not converge unless the coefficient matrix A satisfies the diagonal advantage (in the narrow sense).

2-3-2-2. Gauss-Seidel method implementation

There is no particular change in the graph.

def gauss_seidel(a_matrix, b_array, target_residual):
    """
    Gauss-Seidel method
    Ax = b
    A = (D+L+U)Then
    x^{(m+1)} = D^{-1}(b- L x^{(m+1)} - U x^{(m)})
    x^{(m+1)} = (D+L)^{-1}(b - U x^{(m)})
However, D is a diagonal matrix, L is a lower triangular matrix, and U is an upper triangular matrix.
    
    Parameters
    ----------
    a_matrix : numpy.float64
n × m matrix
    b_array : numpy.float64
n-row matrix
    target_residual : numpy.float64
Positive number. Target residuals.

    Returns
    -------
    x : numpy.float64
n-row matrix
    """
    x_old = b_array.copy()
    lower_matrix = np.tril(a_matrix) #Lower triangular matrix
    upper_matrix = a_matrix - lower_matrix  #Upper triangular matrix
    inv_lower_matrix = np.linalg.inv(lower_matrix)
    count = 0
    while True:
        count += 1
        x = np.dot(inv_lower_matrix, (b_array - np.dot(upper_matrix, x_old)))
        residual = np.linalg.norm(x - x_old) / np.linalg.norm(x)
        x_old = x
        if residual <= target_residual:
            break
        elif count >= 10000:
            import sys
            print(residual)
            sys.exit()
    return x

temperature_gs = temperature_array.copy()
for n in range(Time_step):
    a_matrix = np.identity(len(temperature_gs)) * 2 *(1/d+1) \
                - np.eye(len(temperature_gs), k=1) \
                - np.eye(len(temperature_gs), k=-1)
    temp_temperature_array = np.append(np.append(
                        temperature_lower_boundary, 
                        temperature_gs), temperature_upper_boundary)
    b_array = 2 * (1/d - 1) * temperature_gs + temp_temperature_array[2:] + temp_temperature_array[:-2]
    b_array[0] += temperature_lower_boundary
    b_array[-1] += temperature_upper_boundary
    temperature_gs = gauss_seidel(a_matrix, b_array, 1e-8)
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.plot(x_array, temperature_gs, label="Implicit Gauss-Seidel(100 steps)")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

2-3-3. SOR method

Gauss-Seidel method with relaxation factor $ \ omega $ added. Same as Gauss-Seidel method when $ \ omega = 1 $. When $ A = D + L + U $ using the matrix D with the diagonal components of the matrix A, the upper matrix U, and the lower matrix L,

　Ax = b\\
 (D+L+U) x = b\\
 \tilde{x}^{(m+1)} = D^{-1}(b- L x^{(m+1)} - U x^{(m)}) \quad : Gauss-Seidel\\
 x^{(m+1)} = x^{(m)} + \omega \left(\tilde{x}^{(m+1)} - x^{(m)}  \right)

$ Ax = b $ can be transformed like this, and the method of iteratively finding an approximate solution with the formula at the bottom is called the SOR method (sequential acceleration relaxation method).

x_i^{(m+1)} = x_i^{(m)} + \omega \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(m+1)} - \sum_{j=i}^{n} a_{ij} x_j^{(m)} \right)

2-3-3-1. Implementation of SOR method

This time, the relaxation coefficient $ \ omega $ is a constant. The relaxation factor $ \ omega $ is also available [when the optimal value can be determined](https://ja.wikipedia.org/wiki/SOR method). There is no change in the graph.

def sor(a_matrix, b_array, target_residual):
    """
SOR method
    Ax = b
    A = (D+L+U)Then
    x~^{(m+1)} = D^{-1}(b- L x^{(m+1)} - U x^{(m)}) : Gauss-Seidel
    x^{(m+1)} = x^{(m)} + omega (x~^{(m+1)} - x^{(m)})
However, D is a diagonal matrix, L is a lower triangular matrix, and U is an upper triangular matrix.
    
    Parameters
    ----------
    a_matrix : numpy.float64
n × m matrix
    b_array : numpy.float64
n-row matrix
    target_residual : numpy.float64
Positive number. Target residuals.

    Returns
    -------
    x : numpy.float64
n-row matrix
    """
    x_old = b_array.copy()
    lower_matrix = np.tril(a_matrix) #Lower triangular matrix
    upper_matrix = a_matrix - lower_matrix  #Upper triangular matrix
    inv_lower_matrix = np.linalg.inv(lower_matrix)
    omega = 1.9  #In this example, convergence may be slow.
    count = 0
    while True:
        count += 1
        x_tilde = np.dot(inv_lower_matrix, (b_array - np.dot(upper_matrix, x_old)))
        x = x_old + omega * (x_tilde - x_old)
        residual = np.linalg.norm(x - x_old) / np.linalg.norm(x)
        x_old = x
        if residual <= target_residual:
            break
        elif count >= 10000:
            import sys
            print(residual)
            sys.exit()
    return x

temperature_sor = temperature_array.copy()
for n in range(Time_step):
    a_matrix = np.identity(len(temperature_sor)) * 2 *(1/d+1) \
                - np.eye(len(temperature_sor), k=1) \
                - np.eye(len(temperature_sor), k=-1)
    temp_temperature_array = np.append(np.append(
                        temperature_lower_boundary, 
                        temperature_sor), temperature_upper_boundary)
    b_array = 2 * (1/d - 1) * temperature_sor + temp_temperature_array[2:] + temp_temperature_array[:-2]
    b_array[0] += temperature_lower_boundary
    b_array[-1] += temperature_upper_boundary
    temperature_sor = sor(a_matrix, b_array, 1e-8)
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.plot(x_array, temperature_sor, label="Implicit SOR(100 steps)")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

Summary

The one-dimensional diffusion equation was solved by the central difference method of the explicit method and the implicit method of Crank Nicholson.
In the explicit method of the one-dimensional diffusion equation, it was confirmed that there are restrictions on the lattice width and time step width depending on the diffusion number.
Iterative method was used to implement Crank Nicholson's implicit method.
There are two types of iterative methods: stationary iterative method and unsteady iterative method (Krylov subspace method).
In this implementation, the Jacobi method, Gauss-Seidel method, and SOR method were used.
It was confirmed that the iterative method is not affected by the limitation by the diffusion number as in the explicit method.
The multigrid method and non-stationary iterative method of the stationary iterative method will be explained in the next and subsequent articles.

~~ Next time, I would like to summarize the unsteady iterative method (Krylov subspace method). ~~ I decided to summarize the non-stationary iterative method someday. How to use the unsteady iterative method in Python is summarized in [Python] Article that enables sparse matrix calculation at high speed. Next time, the (linear) convection-diffusion equation, which is a combination of the advection equation of previous article and the diffusion equation dealt with in this article, and its non-linearization I would like to deal with the Burgers equation.

References

https://qiita.com/sci_Haru/items/960687f13962d63b64a0
https://qiita.com/IshitaTakeshi/items/cf106c139660ef138185
http://nkl.cc.u-tokyo.ac.jp/14n/PDE.pdf
https://home.hiroshima-u.ac.jp/nakakuki/other/simusemi-20070907.pdf
https://www.sit.ac.jp/user/konishi/JPN/L_Support/SupportPDF/InplicitMethod.pdf ← Insanely easy to understand
http://ri2t.kyushu-u.ac.jp/~watanabe/RESERCH/MANUSCRIPT/KOHO/GEPP/GEPP.pdf ← Comprehensive and detailed

Iterative reference

The iterative method is the materials and videos listed by Professor Nakajima of the University of Tokyo, and this [site] that summarizes numerical analysis. If you look at (http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?Numerical%20Calculation), you can see it. Professor Nakajima is a god.

Jacobi method references

https://mmm-ssss.com/2019/05/30/jacobi_method_gauss_seidel_method/
https://org-technology.com/posts/solving-linear-equations-yacobi.html
http://www.math.ritsumei.ac.jp/yasutomi/jugyo/Numerical_Analysis/note5.pdf

Steady-state relaxation method implementation references

http://nkl.cc.u-tokyo.ac.jp/13n/SolverIterative.pdf
https://ja.wikipedia.org/wiki/SOR method
The following three websites summarize the implementation of the steady-state relaxation method in a very easy-to-understand manner.
https://org-technology.com/posts/solving-linear-equations-yacobi.html
https://org-technology.com/posts/solving-linear-equations-gauss-seidel.html
https://org-technology.com/posts/solving-linear-equations-sor.html