🎯要点
🎯平流扩散简单离散微分算子 | 🎯相场模拟：简单旋节线分解、枝晶凝固的 | 🎯求解二维波动方程，离散化时间导数
🎯英伟达 A100 人工智能核性能评估模型 | 🎯热涨落流体动力学求解及算法
📜有限差分法 | 本文 - 用例
📜Python和Julia河流湖泊沿海水域特征数值算法模型
📜Python和C++数学物理计算分形热力学静电学和波动方程
📜Python物理学有限差分微分求解器和动画波形传播
📜Python数值和符号算法计算及3D视图物理数学波形方程
📜Python氮氧甲烷乙烷乙烯丙烯气体和固体热力学模型计算
📜Python射频电磁肿瘤热疗数学模型和电磁爆炸性变化统计推理模型
🍇Python有限差分逼近余弦导数
函数
f
(
x
)
f(x)
f(x)在
x
=
a
x=a
x=a点的导数
f
′
(
x
)
f^{\prime}(x)
f′(x)定义为：
f
′
(
a
)
=
lim
x
→
a
f
(
x
)
−
f
(
a
)
x
−
a
f^{\prime}(a)=\lim _{x \rightarrow a} \frac{f(x)-f(a)}{x-a}
f′(a)=x→alimx−af(x)−f(a)
x
=
a
x=a
x=a 处的导数就是此时的斜率。在该斜率的有限差分近似中，我们可以使用点
x
=
a
x=a
x=a 附近的函数值来实现目标。不同的应用中使用了多种有限差分公式，下面介绍其中的三种，其中导数是使用两点的值计算的。
前向差分是使用连接
(
x
j
,
f
(
x
j
)
)
\left(x_j, f\left(x_j\right)\right)
(xj,f(xj)) 和
(
x
j
+
1
,
f
(
x
j
+
1
)
)
\left(x_{j+1}, f\left(x_{j+1}\right)\right)
(xj+1,f(xj+1))的线来估计
x
j
x_j
xj 处函数的斜率：
f
′
(
x
j
)
=
f
(
x
j
+
1
)
−
f
(
x
j
)
x
j
+
1
−
x
j
f^{\prime}\left(x_j\right)=\frac{f\left(x_{j+1}\right)-f\left(x_j\right)}{x_{j+1}-x_j}
f′(xj)=xj+1−xjf(xj+1)−f(xj)
后向差分是使用连接
(
x
j
−
1
,
f
(
x
j
−
1
)
)
\left(x_{j-1}, f\left(x_{j-1}\right)\right)
(xj−1,f(xj−1)) 和
(
x
j
,
f
(
x
j
)
)
\left(x_j, f\left(x_j\right)\right)
(xj,f(xj)) 的线来估计
x
j
x_j
xj 处函数的斜率：
f
′
(
x
j
)
=
f
(
x
j
)
−
f
(
x
j
−
1
)
x
j
−
x
j
−
1
f^{\prime}\left(x_j\right)=\frac{f\left(x_j\right)-f\left(x_{j-1}\right)}{x_j-x_{j-1}}
f′(xj)=xj−xj−1f(xj)−f(xj−1)
中间差分是使用连接
(
x
j
−
1
,
f
(
x
j
−
1
)
)
\left(x_{j-1}, f\left(x_{j-1}\right)\right)
(xj−1,f(xj−1)) 和
(
x
j
+
1
,
f
(
x
j
+
1
)
)
\left(x_{j+1}, f\left(x_{j+1}\right)\right)
(xj+1,f(xj+1)) 的线来估计
x
j
x_j
xj 处函数的斜率：
f
′
(
x
j
)
=
f
(
x
j
+
1
)
−
f
(
x
j
−
1
)
x
j
+
1
−
x
j
−
1
f^{\prime}\left(x_j\right)=\frac{f\left(x_{j+1}\right)-f\left(x_{j-1}\right)}{x_{j+1}-x_{j-1}}
f′(xj)=xj+1−xj−1f(xj+1)−f(xj−1)
为了导出
f
f
f 导数的近似值 ，我们回到泰勒级数。对于任意函数
f
(
x
)
f(x)
f(x)，对于任意函数
f
(
x
)
f(x)
f(x) ，
f
f
f 围绕
a
=
x
j
a=x_j
a=xj 的泰勒级数是
f
(
x
)
=
f
(
x
j
)
(
x
−
x
j
)
0
0
!
+
f
′
(
x
j
)
(
x
−
x
j
)
1
1
!
+
f
′
′
(
x
j
)
(
x
−
x
j
)
2
2
!
+
f
′
′
′
(
x
j
)
(
x
−
x
j
)
3
3
!
+
⋯
f(x)=\frac{f\left(x_j\right)\left(x-x_j\right)^0}{0!}+\frac{f^{\prime}\left(x_j\right)\left(x-x_j\right)^1}{1!}+\frac{f^{\prime \prime}\left(x_j\right)\left(x-x_j\right)^2}{2!}+\frac{f^{\prime \prime \prime}\left(x_j\right)\left(x-x_j\right)^3}{3!}+\cdots
f(x)=0!f(xj)(x−xj)0+1!f′(xj)(x−xj)1+2!f′′(xj)(x−xj)2+3!f′′′(xj)(x−xj)3+⋯
如果
x
x
x 位于间距为
h
h
h 的点网格上，我们可以计算
x
=
x
j
+
1
x=x_{j+1}
x=xj+1 处的泰勒级数以获得
f
(
x
j
+
1
)
=
f
(
x
j
)
(
x
j
+
1
−
x
j
)
0
0
!
+
f
′
(
x
j
)
(
x
j
+
1
−
x
j
)
1
1
!
+
f
′
′
(
x
j
)
(
x
j
+
1
−
x
j
)
2
2
!
+
f
′
′
′
(
x
j
)
(
x
j
+
1
−
x
j
)
3
3
!
+
⋯
f\left(x_{j+1}\right)=\frac{f\left(x_j\right)\left(x_{j+1}-x_j\right)^0}{0!}+\frac{f^{\prime}\left(x_j\right)\left(x_{j+1}-x_j\right)^1}{1!}+\frac{f^{\prime \prime}\left(x_j\right)\left(x_{j+1}-x_j\right)^2}{2!}+\frac{f^{\prime \prime \prime}\left(x_j\right)\left(x_{j+1}-x_j\right)^3}{3!}+\cdots
f(xj+1)=0!f(xj)(xj+1−xj)0+1!f′(xj)(xj+1−xj)1+2!f′′(xj)(xj+1−xj)2+3!f′′′(xj)(xj+1−xj)3+⋯
代入
h
=
x
j
+
1
−
x
j
h=x_{j+1}-x_j
h=xj+1−xj 并求解
f
′
(
x
j
)
f^{\prime}\left(x_j\right)
f′(xj) 得出方程
f
′
(
x
j
)
=
f
(
x
j
+
1
)
−
f
(
x
j
)
h
+
(
−
f
′
′
(
x
j
)
h
2
!
−
f
′
′
′
(
x
j
)
h
2
3
!
−
⋯
)
f^{\prime}\left(x_j\right)=\frac{f\left(x_{j+1}\right)-f\left(x_j\right)}{h}+\left(-\frac{f^{\prime \prime}\left(x_j\right) h}{2!}-\frac{f^{\prime \prime \prime}\left(x_j\right) h^2}{3!}-\cdots\right)
f′(xj)=hf(xj+1)−f(xj)+(−2!f′′(xj)h−3!f′′′(xj)h2−⋯)
括号中的项
−
f
′
′
(
x
j
)
h
2
!
−
f
′
′
′
(
x
j
)
h
2
3
!
−
⋯
-\frac{f^{\prime \prime}\left(x_j\right) h}{2!}-\frac{f^{\prime \prime \prime}\left( x_j\right) h^2}{3!}-\cdots
−2!f′′(xj)h−3!f′′′(xj)h2−⋯ 被称为
h
h
h 的高阶项。高阶项可以重写为
−
f
′
′
(
x
j
)
h
2
!
−
f
′
′
′
(
x
j
)
h
2
3
!
−
⋯
=
h
(
α
+
ϵ
(
h
)
)
-\frac{f^{\prime \prime}\left(x_j\right) h}{2!}-\frac{f^{\prime \prime \prime}\left(x_j\right) h^2}{3!}-\cdots=h(\alpha+\epsilon(h))
−2!f′′(xj)h−3!f′′′(xj)h2−⋯=h(α+ϵ(h))
其中
α
\alpha
α 是某个常数，
ϵ
(
h
)
\epsilon(h)
ϵ(h) 是
h
h
h 的函数，当
h
h
h 变为 0 时，该函数也变为 0。你可以用一些代数来验证这是真的。我们使用缩写“
O
(
h
)
O(h)
O(h)”来表示
h
(
α
+
ϵ
(
h
)
)
h(\alpha+\epsilon(h))
h(α+ϵ(h))，并且一般来说，我们使用缩写“
O
(
h
p
)
O\left(h^p\right)
O(hp)”来表示
h
p
(
α
+
ϵ
(
h
)
)
h^p(\alpha+\epsilon(h))
hp(α+ϵ(h))
将
O
(
h
)
O(h)
O(h) 代入前面的方程得出
f
′
(
x
j
)
=
f
(
x
j
+
1
)
−
f
(
x
j
)
h
+
O
(
h
)
f^{\prime}\left(x_j\right)=\frac{f\left(x_{j+1}\right)-f\left(x_j\right)}{h}+O(h)
f′(xj)=hf(xj+1)−f(xj)+O(h)
这给出了近似导数的前向差分公式为
f
′
(
x
j
)
≈
f
(
x
j
+
1
)
−
f
(
x
j
)
h
f^{\prime}\left(x_j\right) \approx \frac{f\left(x_{j+1}\right)-f\left(x_j\right)}{h}
f′(xj)≈hf(xj+1)−f(xj)
我们说这个公式是
O
(
h
)
O(h)
O(h)。
💦示例：余弦函数 f ( x ) = cos ( x ) f(x)=\cos (x) f(x)=cos(x)。我们知道 cos ( x ) \cos(x) cos(x)的导数是 − sin ( x ) -\sin(x) −sin(x)。尽管在实践中我们可能不知道我们求导的基础函数，但我们使用简单的例子来说明上述数值微分方法及其准确性。以下代码以数值方式计算导数。
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-poster')
%matplotlib inline
h = 0.1
x = np.arange(0, 2*np.pi, h)
y = np.cos(x)
forward_diff = np.diff(y)/h
x_diff = x[:-1:]
exact_solution = -np.sin(x_diff)
plt.figure(figsize = (12, 8))
plt.plot(x_diff, forward_diff, '--', \
label = 'Finite difference approximation')
plt.plot(x_diff, exact_solution, \
label = 'Exact solution')
plt.legend()
plt.show()
max_error = max(abs(exact_solution - forward_diff))
print(max_error)
0.049984407218554114
如上图所示，两条曲线之间存在微小的偏移，这是由于数值导数求值时的数值误差造成的。两个数值结果之间的最大误差约为 0.05，并且预计会随着步长的增大而减小。
💦示例：以下代码使用递减步长 h h h 的前向差分公式计算 f ( x ) = cos ( x ) f(x)=\cos (x) f(x)=cos(x) 的数值导数。然后，它绘制近似导数和真实导数之间的最大误差与 h h h 的关系，如生成的图形所示。
h = 1
iterations = 20
step_size = []
max_error = []
for i in range(iterations):
h /= 2
step_size.append(h)
x = np.arange(0, 2 * np.pi, h)
y = np.cos(x)
forward_diff = np.diff(y)/h
x_diff = x[:-1]
exact_solution = -np.sin(x_diff)
max_error.append(\
max(abs(exact_solution - forward_diff)))
plt.figure(figsize = (12, 8))
plt.loglog(step_size, max_error, 'v')
plt.show()
双对数空间中直线的斜率为 1 ；因此，误差与 h 1 h^1 h1成正比，这意味着，正如预期的那样，前向差分公式为 O ( h ) O(h) O(h)。
