🎯要点
🎯平流扩散简单离散微分算子 | 🎯相场模拟：简单旋节线分解、枝晶凝固的 | 🎯求解二维波动方程，离散化时间导数
🎯英伟达 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)。