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

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

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

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

---

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

用向后差分公式近似Ut，

用中心差分公式近似Uxx，

用差分公式替代在点的热方程，得到

其中令 

则上述方程可以写成

改写成的矩阵方程

其中

则

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

---

偏微分方程编程步骤

1. 导入numpy，matplotlib包
2. 定义初值函数
3. 设置空间和时间步数，以及定义域范围
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('T', fontsize=14)
ax1.set_ylabel('X', 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()

---

示例运行结果