📜有限差分-用例
📜离散化偏微分方程求解器和模型定型 | 📜三维热传递偏微分方程解 | 📜特定资产期权价值偏微分方程计算 | 📜三维波偏微分方程空间导数计算 | 📜应力-速度公式一阶声波方程模拟二维地震波 | 📜微磁学计算磁化波动求解器、色散关系和能垒的弦法 | 📜磁倾斜导数数据平滑
📜指数衰减：🖊常微分方程数值求解器 | 🖊绘制衰减图 | 🖊绘制（正向欧拉、反向欧拉和克兰克-尼科尔森）西塔规则算法放大因子图 | 🖊泰勒级数展开符号计算三种算法误差 | 🖊模型误差、数据误差、离散化误差和舍入误差 | 🖊求解器泛化
📜Python热涨落流体力学求解算法和英伟达人工智能核评估模型
📜常微分方程用例：Python机器人动力学和细胞酶常微分方程
✒️Python不同初始条件下热方程
有限差分法是获得偏微分和代数方程数值解的技术之一。在该方法中，解在有限网格点中以离散形式近似。
首先考虑一个偏微分方程：
u
t
+
a
u
x
=
0
u_t+a u_x=0
ut+aux=0
正向时间前向空间算法由下式给出：
V
m
n
+
1
−
V
m
n
k
+
a
V
m
+
1
n
−
V
m
n
h
=
0
\frac{V_m^{n+1}-V_m^n}{k}+a \frac{V_{m+1}^n-V_m^n}{h}=0
kVmn+1−Vmn+ahVm+1n−Vmn=0
正向时间中心空间算法由下式给出：
V
m
n
−
1
−
V
m
n
k
+
a
⋅
V
m
−
−
−
V
m
−
1
n
2
h
−
0
\frac{V_m^{n-1}-V_m^n}{k}+a \cdot \frac{V_{m-}^{-}-V_{m-1}^n}{2 h}-0
kVmn−1−Vmn+a⋅2hVm−−−Vm−1n−0
中心时间中心空间算法由下式给出
V
m
n
+
1
−
V
m
n
−
1
2
k
+
a
⋅
V
m
+
1
n
−
V
m
−
1
n
2
h
=
0
\frac{V_m^{n+1}-V_m^{n-1}}{2 k}+a \cdot \frac{V_{m+1}^n-V_{m-1}^n}{2 h}=0
2kVmn+1−Vmn−1+a⋅2hVm+1n−Vm−1n=0
让我们考虑另一个偏微分方程，
u
t
=
b
u
x
x
;
b
>
0
u_t=b u_{x x} ; \quad b>0
ut=buxx;b>0
正向时间中心空间算法由下式给出：
V
m
n
+
1
−
V
m
n
k
=
b
V
m
+
1
n
−
2
V
m
n
+
V
m
−
1
n
h
2
=
0
\frac{V_m^{n+1}-V_m^n}{k}=b \frac{V_{m+1}^n-2 V_m^n+V_{m-1}^n}{h^2}=0
kVmn+1−Vmn=bh2Vm+1n−2Vmn+Vm−1n=0
示例：数值求解
u
t
=
0.05
u
x
x
u_t=0.05 u_{x x}
ut=0.05uxx
- u u u 代表温度
- x x x 表示 0 ≤ x ≤ L 0 \leq x \leq L 0≤x≤L 的位置
- t t t 表示 t > 0 t>0 t>0的时间
- 边界条件为 u ( t , 0 ) = 0 u(t, 0)=0 u(t,0)=0 和 u ( t , L ) = 0 u(t, L)=0 u(t,L)=0（ t > 0 t>0 t>0）
- 初始条件为 u ( 0 , x ) = sin ( π x ) u(0, x)=\sin (\pi x) u(0,x)=sin(πx) 对于 0 ≤ x ≤ L 0 \leq x \leq L 0≤x≤L
- b b b 表示 b > 0 b>0 b>0 的扩散系数
代码求解：
import numpy as np
import matplotlib.pyplot as plt
L = 1
T = 1
m = 5
n = 5
h = L / m
k = T / n
b = 0.05
mu = k / h**2
c = b * mu
if c <= 0 or c >= 0.5:
print('Scheme is unstable')
v = np.zeros((m + 1, n + 1))
ic1 = lambda x: np.sin(np.pi * x)
for j in range(1, m + 2):
v[0, j - 1] = ic1((j - 1) * h)
b1 = lambda t: 0 # L.B.C
b2 = lambda t: 0 # R.B.C
for i in range(1, n + 2):
v[i - 1, 0] = b1((i - 1) * k)
v[i - 1, n] = b2((i - 1) * k)
for i in range(n):
for j in range(1, m):
v[i + 1, j] = (1 - 2 * b * mu) * v[i, j] + b * mu * v[i, j + 1] + b * mu * v[i, j - 1]
x = np.linspace(0, L, m + 1)
t = np.linspace(0, T, n + 1)
X, T = np.meshgrid(x, t)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, T, v, cmap='viridis')
ax.set_xlabel('Space X')
ax.set_ylabel('Time T')
ax.set_zlabel('V')
plt.title('Python for Heat')
plt.show()
接上例，初始条件改为：
对于 0 ≤ x ≤ L 0 \leq x \leq L 0≤x≤L
u ( 0 , x ) = { 2 x if x < 0.5 2 ( 1 − x ) 否则 u(0, x)= \begin{cases}2 x & \text { if } x<0.5 \\ 2(1-x) & \text { 否则 }\end{cases} u(0,x)={2x2(1−x) if x<0.5 否则
代码数值解：
import numpy as np
import matplotlib.pyplot as plt
L = 1
T = 1
m = 5
n = 5
h = L / m
k = T / n
b = 0.05
mu = k / h ** 2
c = b * mu
if c <= 0 or c >= 0.5:
print('Scheme is unstable')
v = np.zeros((m + 1, n + 1))
ic1 = lambda x: 2 * x
ic2 = lambda x: 2 * (1 - x)
x = np.linspace(0, L, m + 1)
x = np.linspace(0, L, m + 1)
for j in range(1, m + 2):
if x[j - 1] < 0.5:
v[0, j - 1] = ic1(x[j - 1])
else:
v[0, j - 1] = ic2(x[j - 1])
b1 = lambda t: 0 # L.B.C
b2 = lambda t: 0 # R.B.C
for i in range(1, n + 2):
v[i - 1, 0] = b1((i - 1) * k)
v[i - 1, n] = b2((i - 1) * k)
for i in range(n):
for j in range(1, m):
v[i + 1, j] = (1 - 2 * b * mu) * v[i, j] + b * mu * v[i, j + 1] + b * mu * v[i, j - 1]
# Visualization
x = np.linspace(0, L, m + 1)
t = np.linspace(0, T, n + 1)
X, T = np.meshgrid(x, t)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, T, v, cmap='viridis')
ax.set_xlabel('Space X')
ax.set_ylabel('Time T')
ax.set_zlabel('V')
plt.title('Python for Heat ')
plt.show()
