Numerical Processing | XiangSiqi’s Blog

type

Post

status

Published

date

Nov 12, 2025

slug

CP4

summary

Interpolation

Lagrange Interpolation

从一个简单的例子出发，考虑点和，那么直线为：

考虑，显然满足，我们称这样的组合为线性插值基函数。

我们定义插值函数为两个插值基函数的线性组合，系数为对应点上的函数值。

因此我们有二次插值函数，它的插值基是不超过二次的多项式，满足：

即：

组合出：

对于任意二次多项式，必然可以展开为，由利用得到：

这就是拉格朗日多项式：

这是一个次多项式，对应于个点，插值函数为：

拉格朗日插值的结果有无穷多阶导数，这是激进的。而分段线性函数只有1阶导数，但是在采样点上转折，这也是激进的。

同时，拉格朗日插值法的零阶项在各个点都相等，但是1阶及以上都不相等，而泰勒展开在展开点满足各阶导数相等。

我们需要找到一个折中的思路。

Hermite interpolation

Hermite Polynomial

我们首先直接给出结论，对于函数阶可导，它能近似为厄米多项式，即：

这里:

其中为采样的点，而：

显然，对于零阶项，当时，，，因此：

而对于1阶项，可以证明当时：

而对于第二项：

这就证明了一阶导数：

厄米多项式方法意味着，对于个样本点，逼近函数就是阶可导的。对于任意的个样本点，，都可以找到一个三阶多项式去逼近它。

Pchip

实际中我们不使用这个方法，而是采用所谓的分段三次Hermite插值，它的思路是，分段考虑的两个样本点，取。考虑构造多项式：

考虑条件为，则有：

Spline

除了这个方法外，还可以控制一阶和二阶导数都相等，这就是spline方法（三次样条插值函数）。

考虑在子区间是，考虑零阶条件：

总共贡献个方程，考虑连续条件：

其实这一条被上一条包含，考虑一阶导数：

有个方程，二阶导数：

有个方程，因此以上总共个方程，但是总共有个未知数，剩下两个方程来源于边界条件，通常采取自然边界：

或者固定边界：

或许这种方法更加实用，因为它不要求提供导数信息。例子如下

Scipy Interpolate

首先是样条

kind可以选择linear对应于线性插值，zero,slinear,quadratic,cubic对应于0，1，2，3阶样条插值，nearest对应于直接选择要插值的点的最近点,如下图

三阶样条如下：

然后是厄米插值：

这个函数直接返回插值所对应的值，也有返回函数的：

Curve Fit

Least Squares Method

考虑插值结果，真实值为，偏差为，有一下三种拟合问题：

使得绝对值和最小：

使得偏差绝对值的最大值最小

使得偏差的平方和最小(LSM)

或者加权平方和：

Polynomial Fit

考虑插值函数为：

从而将成本函数写为：

考虑求导数：

从而为：

可以容易看到这是一个矩阵方程，

可以直接将系数解出来

Polynomial Fit Command

这里deg是多项式的阶数，函数返回的是多项式的系数。

从高次到低次。而要调用系数：

或者获得函数对象

Exponential Fit

指数函数形式为：

取对数：

这样就变成了一个线性拟合问题，只需要将数据的取对数，然后和进行线性拟合（1阶最小二乘法），得到的截距取指数，即可得到。

General Form of Linear Least Squares

取在上线性无关的基函数，拟合函数为：

偏差为：

成本函数为：

考虑求导数：

可以得到：

其中，这样就仍然可以写成矩阵方程式：

我们选取正交基，这样系数矩阵为对角阵，这样：

这样最小二乘函数：

Scipy Curve Fit

基本函数是

这里popt是拟合参数，pcov是参数的协方差矩阵

Linear Fit

结果是

Exponential Fit

结果如下：

Gauss Function Fit

Fourier Series Fit

结果为：

Zeros of The Equation

Monotonic function

Divide Method

对分法的思路是，如果找到单调函数的，那么零点一定在之间，只需要考虑是否为正，若为正，则零点在之间，若为负则在之间，从而继续迭代即可。我们用这个方法来解分子振动的半经典量子化。

Tangent Method

考虑在处泰勒展开：

并且考虑时，展开值为0，即，则：

这表明是点的切线和轴的交点。

我们的思路就是反复计算这个点，从而逼近零点。

但是这种方法要求函数必须可导，为了解决这个问题，我们考虑弦割法。

Wire Cut Method

我们把迭代公式的导数换成弦的割线斜率：

示意图如下：

不断迭代即可。

Scipy, Numpy, Sympy

scipy.optimize.fsolve

这里的数字0是初始值，对于多根的方程，从不同的初始值出发得到不同的根

scipy.optimize.brentq

采用割弦法求解范围。

这里传入的两个数字是割弦法的两个端点，因此根只会在两个端点中

numpy.poly1d

poly1d返回的对象有两个属性.r返回的是多项式的根，.c返回的是多项式的系数

sympy解方程组

fsolve解方程组

值得注意的是，传入的三个仍然是三个变量的初值

Fast Fourier Transform

Discrete Fourier Transform

标准傅立叶变换为：

该傅立叶变换提取中的成分。考虑为对应周期，则测量时间中能够测量到个波，考虑离散化为，这里为总测量次数，得到：

这就是离散形式的傅里叶变换，取，得到：

对逆变换采取同样的处理得到：

考虑取得到：

这里.

显然离散形式变换可以写成矩阵形式：

指数项，这里代表测量时间内波的数目，因此意义上是频率，意义上是圆频率，尽管实际上我们称为时间，为时间意义上的频率。

逆变换强度项相当于各频率的振幅，振幅的平方称为功率谱

，即零频分量是时域分量之和，为尼奎斯特频率，为采样频率的一半。除了零频和尼奎斯特频率，其他频率两两共轭。

Fourier Series and Integration

标准的傅里叶级数，对于定义在的任意平方可积函数可以展开为：

有复数指数形式：

而傅立叶积分为：

Numpy FFT

来自numpy.fft模块

一维FFT及其逆变换：

得到的结果就是变换得到的一维数组

二维FFT及其逆变换

得到的结果就是变换得到的二维数组

n维FFT及其逆变换

得到的结果就是变换得到的n维数组

得到数组后，我们还需要频率绘制出，需要：

离散傅立叶变换频率

即返回. 因此所有FFT数据都是按照从0,1,2,3,..,(N-1/2),-(N-1/2),..,-3,-2,-1（奇数），0,1,2..,(N/2),-(N/2-1),-2,-1（偶数）