用Python计算偏微分方程数值解——隐式方法

偏微分方程主要有三类：椭圆方程，抛物方程和双曲方程。

本文采用隐式有限差分法求解偏微分方程，通过典型案例（热方程）来演示隐式方法的使用。

相比一般的前向差分方法，隐式方法对于所有选择的步长都无条件稳定。

---

一维热传导方程的数学模型

$\begin{cases} u_t = Du_{xx}\\ u(x,0) = sin^22\pi x &\text{初值条件}\\ u(0,t)=0 &\text{边值条件}\\ u(1,t)=0 \\ \end{cases}$

用向后差分公式近似$u_t$，

$u_t = \frac{1}{k}(u(x,t)-u(x,t-k)) + \frac{k}{2}u_{tt}(x,c_0)$

用中心差分公式近似$u_{xx}$，

$u_{xx}=\frac{1}{h^2}(u(x+h,t)-2u(x,t)+u(x-h,t))$

用差分公式替代在点$(x_i,t_i)$的热方程，得到

$\frac{1}{k}(\omega_{i,j} - \omega_{i,j-1}) = \frac{D}{h^2}(\omega_{i+1,j} - 2\omega_{i,j} + \omega_{i-1,j})$

其中令$\sigma = \frac{Dk}{h^2}$

则上述方程可以写成

$-\sigma\omega_{i+1,j} + (1+2\sigma)\omega_{i,j}-\sigma\omega_{i-1,j} = \omega_{i,j-1}$

改写成$m\times m$的矩阵方程

$\begin{bmatrix} 1+2\sigma & -\sigma & 0 & \cdots & 0 \\ -\sigma & 1+2\sigma & -\sigma & \ddots & 0 \\ 0 & -\sigma & 1+2\sigma & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & -\sigma & 1+2\sigma\\ \end{bmatrix} \begin{bmatrix} \omega_{1,j} \\ \vdots \\ \omega_{m,j} \\ \end{bmatrix} = \begin{bmatrix} \omega_{1,j-1} \\ \vdots \\ \omega_{m,j-1} \\ \end{bmatrix} + \sigma \begin{bmatrix} \omega_{0,j} \\ 0 \\ \vdots \\ 0\\ \omega_{m+1,j} \\ \end{bmatrix}$

其中

$A = \begin{bmatrix} 1+2\sigma & -\sigma & 0 & \cdots & 0 \\ -\sigma & 1+2\sigma & -\sigma & \ddots & 0 \\ 0 & -\sigma & 1+2\sigma & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & -\sigma & 1+2\sigma\\ \end{bmatrix}$

则$\omega_j = A^{-1}\omega_{j-1} + b$

只需用$j-1$时刻的数值不断迭代，就可以获得所有点上的数值解

---

偏微分方程编程步骤

1. 导入numpy，matplotlib包
2. 定义初值函数
3. 设置空间和时间步数，以及定义域范围$(0,1)$
4. 设置矩阵A，递推求解方程在区间内的数值解
5. 画图

---

Python示例

# -*- coding: utf-8 -*-
"""
Created on Sat Jul  6 08:18:27 2024

@author: chlj
"""

# Solving partial differential equations
# 偏微分方程数值解法

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# 抛物线方程
# U_t = D * U_xx
# 隐式方法边界条件U(0,t) = 0, U(1,t) = 0

# 初始条件函数 U(x,0)
def funcUx0(x): 
    u0 = (np.sin(2 * np.pi * x)) ** 2
    return u0

M = 101  #空间步数
N = 101  #时间步数
h = 1/(M-1)  #x的步长
k = 1/(N-1)  #t的步长

D = 1  #扩散系数
sigma = D*k /(h*h)  #σ

#设置矩阵
A = np.zeros((M-2, M-2))
w = np.zeros((M-2, N))
for i in range(0,M-2):
    for j in range(0,M-2):
        if j == i:
            A[i][j] = 1 + 2 * sigma
        elif j == i+1:
            A[i][j] = -sigma
        elif j == i-1:
            A[i][j] = -sigma
'''
        [[ 201. -100.    0. ...    0.    0.    0.]
         [-100.  201. -100. ...    0.    0.    0.]
         [   0. -100.  201. ...    0.    0.    0.]
         ...
         [   0.    0.    0. ...  201. -100.    0.]
         [   0.    0.    0. ... -100.  201. -100.]
         [   0.    0.    0. ...    0. -100.  201.]]
'''

A_inv = np.linalg.inv(A) #求逆

b = np.zeros((M-2, 1))
for i in range(1,M-1):
    b[i-1][0] = funcUx0(k*i)
w[:,0]= b.T

for i in range(0,N-2):
    b = np.dot(A_inv,b)
    w[:,i+1] = b.T
    
c = np.zeros((1,N)) #边界条件

w = np.r_[c,w,c]  #添加边界条件

#画网格
x, t = np.meshgrid(np.linspace(0, 1, M),
				   np.linspace(0, 1, N))

#画图
f = plt.figure('Wireframe', facecolor='lightgray')
ax1 = f.add_subplot(111,projection='3d')
ax1.set_xlabel('X', fontsize=14)
ax1.set_ylabel('T', fontsize=14)
ax1.set_zlabel('Z', fontsize=14)
ax1.set_zlim(-1, 1)
ax1.plot_wireframe(x, t, w, linewidth=0.5)
ax1.view_init(elev=20, azim=-120)
plt.savefig('抛物线方程.jpg',dpi = 200)
plt.show()

示例运行结果