Numerical Differentiation and Numerical Integration

type

Post

status

Published

date

Nov 11, 2025

slug

CP3

summary

Numerical Differentiation and Numerical Integration

Numerical Differentiation

Central Difference

向前和向后差分泰勒展开：

得到：

我们对向前和向后差分取平均，得到中心差分：

考虑误差，将泰勒展开到更高阶：

相减：

考虑二阶导数：

这是向前差分，代入向后差分的一阶导数：

考虑误差评估，泰勒展开到4阶：

从而：

我们得到如下：

Lagrange Polynomial

考虑对原函数采样的点为，考虑找到来近似原函数，并且，利用拉格朗日多项式可以做到这件事：

显然，拉格朗日多项式满足，这意味着对于求和：

这就满足了我们的要求。但是我们的问题是，这样一个表达式和原函数究竟差多少呢？我们先给出结论：

现在我们来证明，考虑构造：

由于拉格朗日多项式存在，这个函数必然满足：

并且显然有导数

又是阶多项式，导数为0，考虑最后一项，分母为常数，分子为的阶多项式，且系数为1，因此结果为，得到：

根据罗尔中值定理，必定存在，这里，使得：

得到：

Three Point Formular

一阶导数

现在考虑只有三个点，且，则根据上一节的公式：

根据定义：

求导数：

分别代入：

这就是三点一阶导数公式。

二阶导数

再求导数：

注意到0阶项：

就是中点公式。

Rounding Error and Discrete Error

考虑中心差分的一阶导数公式：

这里的都是浮点数计算，有一个舍入误差：

从而：

完整的误差项：

第一项称为舍入误差，第二项称为离散误差，我们设,则：

对其求导数：

这表明在取到极小值，表明在有舍入误差的情况下，采样不是越密越好，这是因为过密会导致数据差值下降到和误差一个量级，数据完全失真。

Numerical Integration

Trapezoid Formula

我们可以容易给出：

这里同样的问题是向前还是向后取值：

这里。考虑取平均：

其实就是中点的乘以间距：

因此可以不等步长采样。

Simpson Formular

考虑三点积分，那么积分：

将函数在展开：

观察得到，奇阶项的积分都为0，因此：

现在考虑用中心二阶差分公式：

代入：

两个四阶项近似消去：

实际计算使用这个表达式迭代求和即可，现在考虑对积分全程求和：

这就是所谓的Simpson Formular.

Numpy and Scipy

numpy diff

findiff

它的安装方法是：

使用方法是：

numpy gradient

还可以进行多方向梯度：

此时会返回两个梯度矩阵

Laplace

容易给出：

Scipy Trapezoid

直接调用：

这里的x是用于计算非均匀间隔，而dx则是用于当没有x时采取的默认间隔，即原公式的

还有一个函数，cumulative_trapezoid用于返回每个点的积分值：

这里initial放在输出矩阵的第一位用于保持矩阵形状不变。

Scipy Simpson

同样直接调用：

根据理论，simpson方法要求传入奇数个点，当传入数列不满足时，由even参数控制行为。当even = arg意味着传入前N-2个点和后N-2个点的结果的平均值，如果是even = first/last则只取其中一个。

可以看到，simpson方法精度更高。

Scipy quad and dblquad

直接调用即可

这里epsabs和epsrel分别为绝对误差和相对误差，计算会在误差小于这两者之后停止。

而对于二元积分，通常写为：

它使用dblquad方法计算：

对于三元积分：

它使用tplquad方法计算：

它们都可以指定误差。

Sympy Symbolic Computation

Symbolic Object

可以使用symbols类：

或者使用sympy.abc模块：

Limit

直接使用极限函数即可：

如果两侧值不同，例如在两侧不同，则返回

Taylor Expansion

直接调用series函数即可：

传入函数，变量x，在x0展开到n阶，返回的是符号表达式：

Differentiation and Integration

直接调用函数即可：

分别是函数，自变量，阶数

积分同样：

如果传入的是三元元组(x,a,b)，则是定积分

Solve Equation

首先创建Equation，然后直接调用solve函数即可：

Examples

Magnetic Field of a Ring Current

首先建模，考虑电流环的电流元的坐标：

这里是电流环的半径，是环上角度。考虑任意点坐标，则根据Biot-Savart定律：

这里：

在代码中实现：

这里用CoordSys3D创建3维空间，并创建向量，进行表达式计算。并利用lambdify转化成函数，在循环中将值传入转化为1元函数并进行积分，最后利用quiver打印。

Semiclassical Quantization of Molecular Vibrations

考虑势能是Lennard-Jones势：

它的图像是：

考虑能量守恒：

考虑，则取到两个半径，左侧称为，右侧称为，在这两处，这意味着相空间的闭合曲线关于x轴对称。

考虑波尔的轨道量子化：

回路积分因此可以切成两半：

考虑化简为：

其中：

同理，和满足方程。因此我们必须迭代各个，解出，并逐步增大的以找到积分成立的。

这里使用sympy直接把的显式依赖解出来加快计算。