差分法求解微分方程的示例：
代码语言：javascript复制import numpy as np
from numpy import exp
from sympy import symbols, solve, Eq
from math import e
from matplotlib import pyplot as plt
# d2phi/dx2 d phi/dx - 2* phi =0, phi(0) = 0, phi(4)=1
# 通用部分开始
x_l, p_x_l = 0, 0 # 左端 边界条件
x_r, p_x_r = 4, 1 # 右端 边界条件
n = 17 # 节点数, n>=3
X = np.linspace(x_l, x_r, n)
dx = X[1] - X[0]
P = [symbols(f'p{i}') for i in range(n)] # 定义 p0, p1, ..., pn # p for phi
Dp_dx = [symbols(f'dp_dx{i}') for i in range(n)] # 一阶导数
D2p_dx2 = [symbols(f'd2p_dx2{i}') for i in range(n)] # 二阶导数
# print(P)
# print(X)
# 左边界处
Dp_dx[0] = (-3 * P[0] 4 * P[1] - P[2]) / (2 * dx)
D2p_dx2[0] = (P[0] - 2 * P[1] P[2]) / (dx * dx)
# 中间节点
for i in range(1, n-1):
Dp_dx[i] = (P[i 1] - P[i-1]) / (2 * dx)
D2p_dx2[i] = (P[i 1] - 2*P[i] P[i-1]) / (dx * dx)
# 右边界
Dp_dx[n-1] = (3*P[n-1] - 4*P[n-2] P[n-3]) / (2 * dx)
D2p_dx2[n-1] = (P[n-1] - 2 * P[n-2] P[n-3]) / (dx * dx)
# 通用部分结束
eqs = [Eq(P[0], p_x_l), Eq(P[n-1], p_x_r)] # 边界条件
for i in range(1, n-1): # 内部满足微分方程
eq = Eq(D2p_dx2[i] Dp_dx[i] - 2 * P[i], 0)
eqs.append(eq)
# print(eqs)
r = solve(eqs,P)
# print(r)
p = np.zeros(n)
for i in range(n):
p[i] = r[P[i]]
print(p)
x_ = np.linspace(x_l, x_r,1000)
y = (np.exp(x_) - np.exp(-2*x_)) / (e**4 - e**(-8)) # 真实解
plt.plot(x_, y, c="g", label="真实解")
x = np.linspace(x_l, x_r, n)
plt.scatter(x, p, c="r", label="差分解")
plt.legend(loc="lower right")
plt.show()
附表：
