cdinten的专栏

前面可以看到，通过使用2个到4个不同的等距节点，采用Lagrange插值公式构造多项式曲线拟合原待求积分的函数，实现数值积分，从两点的梯形公式，到三点的1/3 Simpson公式，以及四点的3/8 Simpson公式，其本质都是一样的，按这种思路，我们可以一直做下去，例如五点、六点等的积分公式，那么，这其中的方法，就是Newton-Cotes公式。由前面介绍的曲线拟合中的插值算法，可知这些点及其对应函数值，通常可以构成次数不超过n次的多项式

既然前的Simpson可以通过使用三个点构造二次曲线近似积分，那么，如果点数增加到了4个，然后不就可以构造三次多项式的曲线，实现对目标值的积分吗？如果采用和上一节介绍的同样的方法，也是可以实现的。

才想起来昨天是1024，程序员也有自己的节日了。

前面我们通过用矩形和梯形的数值算法，近似实现了数值积分，那么，和之前插值类似，是否可以使用多项式来拟合曲线，然后将该多项式作为被积函数求积分呢？当然是可行的，如果以最简单的二次多项式为例，如果给出三个点及对应的函数值，一般就可以构造出一个二次多项式，然后以此进行积分，这就是Simpson（辛普森）积分公式的思想。

梯形积分公式其实和矩形积分公式类似，只是前面在采用矩形积分公式时，使用的是中点处的纵坐标作为矩形的高，而在梯形积分公式中，使用的是边界上两个点的高度，作为直角梯形的上底和下底，然后横坐标的差作为梯形的高，其近视积分为：可以看到，如果对于线性函数，这种算法是精确的。

前面介绍了数值积分的基本背景知识，接下来就介绍各种常见的数值积分算法，本次主要介绍矩形和梯形积分公式。

众所周知，对于一个函数求微分，只要搞清楚了链式法则，对已知解析式的函数，都可以比较方便的求出其导数的解析式。但是反过来，即使是知道一些函数的解析式，也不一定（或者说大概率）没法算出其积分的解析式，这样的函数简直不胜枚举，例如等等，甚至，在有些时候连原函数解析式是怎样的都不知道，只有一系列的采样点，而又要通过这些采样点进行积分计算。通过数值方法计算的基本方法叫做数值积分，通常使用来近似。

理查森外推法( Richardson extrapolation)是一种提高某些数值过程精度的简单方法，在数值方法中广泛应用。理查森外推法的基本思想是通过对原函数进行多次求导，并在每一步求导的基础上进行线性组合，得到一个新的函数，这个新的函数与原函数的差距会逐渐减小，最终趋于零。理查森外推法的优点是简单易行，适用于各种数值计算。理查森外推法的一般步骤如下：1. 确定原函数的导数。2. 确定原函数的导数的导数。3. 利用导数的导数，构造新的函数。

根据前面的定义，虽然我们无法通过数值计算取得极限为零，但是当变化量足够小的时候，是可以通过差分的方法，估算处该点处的导数值的，通过对泰勒多项式的展开，以及差分选用的方式，分为了前向差分，中间差分和后向差分，以下是误差限为。所以，在数值计算的时候，虽然公式可以保证，较小的步长可以产生更加精确的结果，但是由于在实际计算时，由于浮点数计算的限制，以及分母太小，容易造成分子的误差被放大，最终结果是更小的步长并不一定会产生更小的误差。处的导数，就是该点处切线的斜率，导数对于函数的变化率至关重要。

Python虽然作为一种解释性语言，通常是使用解释器处理脚本的方式完成工作，但是也具有和很多其他诸如C、C#语言等类似的方式，可以将.py文件编译为*.pyd文件，即Python动态链接库文件。关于动态链接库的优缺点这里就不展开细说了，具体介绍如何编译及使用。要将Python文件编译成.pyd文件，可以使用Cython，以下是基本的操作步骤。

近日将Python从3.11升级到了3.12，然后把几个包例如numpy等也通过pip给upgrade了一下，结果原来运行的好好的脚本，都运行不了了，还出现各种报错。怀疑是自己升级了环境导致的，因此通过搜索引擎检索了一下，有这么个方式可以解决，做个笔记，也和大家分享一下。注意：保险起见，在升级所有包之前，建议备份你的环境或确保重要项目在一个独立的虚拟环境中，避免更新导致的兼容性问题。既然如此，那就通过命令，再安装一个pip-review的包。然后就可以正常运行之前的脚本了。可以发现很多过时的软件包。

一开始，查找到的资料都是使用imshow，但是这会有带来几个问题，一个是图形的原点发生了变化，另外一个问题就是图形比例也产生了变化，当然最大的问题是图形占据了整个绘图区域，完全喧宾夺主了，与我们设想的只在绘图区域中占据很小的一块不相符。在使用matplotlib下的pyplot绘图时，有时处于各种原因，需要采用类似贴图的方式，插入外部的图片，例如添加自己的logo，或者其他的图形水印等。就是左上方了，其中的四个元素分别表示的是归一化之后，子绘图区域占据父绘图区域的x、y，以及width和height。

PCHIP是“分段三次埃尔米特插值多项式（piecewise cubic Hermite interpolating polynomial）”的缩写，但是单纯根据这个名字，其实并不能确定是哪一种分段三次埃尔米特插值多项式，因为样条插值也是分段三次埃尔米特插值多项式的一种，只是对斜率的约束条件不同而已，这里讨论的是一种保形（shape preserving）的特定插值函数，和Akima插值算法类似，其关键的思想还是如何确定各个采样点的斜率，使函数的值不会再局部上过多的偏离给定的数据。

在前面的Akima插值中，计算斜率使用如下公式：如果记：在出现分母分子同时为零的情况时，会出现NaN的计算结果，Akima他自己也意识到这种问题，因此，在原来的算法上做了修订，也就是Modified Akima算法，有时也被简写为makima算法。修订如下：这样的话，就解决了在出现超过2个点具有连续值时的无法被零除问题，此时。

前面介绍到的很多插值算法，在多项式阶次较高时，都会存在龙格现象，最直观的表现就是导数的变化很急剧，所以一般采用的方式是使用分段的三次插值，例如前面介绍到的三次样条曲线，这里再介绍Akima插值算法，通过该方法生成的插值曲线，具有更加平滑的特性。

前面介绍到紧固和自然三次样条曲线，这次介绍一下非扭结点三次样条曲线。所谓的非扭结点，是指由于最开始的两个子区间使用插值多项式相同，最后两个子区间所使用的插值多项式也相同，这就会导致在这段多项式上起不到扭结点的效果，故而得名not-a-knot。

前面介绍到紧固三次样条曲线，这次介绍一下自然三次样条曲线。

前面介绍到了三种样条曲线的类型，接下来将分三个章节分别介绍紧固三次样条曲线，自然三次样条曲线和非扭结点三次样条曲线。

介绍三次样条插值的基本概念

这次介绍一下使用差商表构造Hermite多项式的方法。前面介绍到了两种很经典的插值多项式，即Lagrange和Newton插值多项式，并在前一篇中阐述了如何通过Lagrange插值方式构造Hermite多项式，这次通过牛顿差商法构造Hermite多项式。

matplotlib之动态拖动点。

在scipy中，scipy.interpolate下还有一个PPoly的类，用于表示插值多项式，很多插值算法的结果，都以该类的实例返回，因此有必要了解该类的使用方法。

不管是前面介绍到拉格朗日插值还是牛顿插值，拟合的函数比线性插值更加“优秀”，即它们都是连续可导的，但是，有时拟合还有这样的要求，就是除了在给定点处的函数值要相等外，还要求在这些指定点处的导数值也要是制定的值，甚至是m阶导数相等，也就是所谓的密切多项式，处理这类拟合的插值方法，就是Hermite插值。

插值是根据已知的数据序列（可以理解为你坐标中一系列离散的点），找到其中的规律，然后根据找到的这个规律，来对其中尚未有数据记录的点进行数值估计的方法。最简单直观的一种插值方式是线性插值，它是数学、计算机图形学等领域广泛使用的一种简单插值方法，根据已知数据点来估某一点的函数值。线性插值法的特点是简单、方便，适用于函数在某一段区间内是线性的情况，即函数在该区间内可以用一条直线来近似表示。线性插值法的优点是计算量小，缺点是对于。多个点而言在连接点出出在明显的间断点。

为了便于自己和大家查阅，将链接都集中在这里，不定时更新。

总结Numpy中数据的保存和载入

虽然在给定了N个点以后，通过这个点的最小幂多项式是确定的，但是表达方式可不止一种，例如前面提到的系数方式，根方式，还有插值的Lagrange形式等。求解差商的秘诀还是构造差商表，和前面Neville方法类似，但是最后的构造形式会有些差别。使用这种表达方式的关键自然是确定前面的系数，而Newton差商公式，就是解决这一问题的。显然这个式子最高次幂不会超过n，然后在。零阶差商用函数值表示，最后所求的系数。显然这和之前的多项式是一致的。在差商矩阵的主对角线上。

上图汇总，红色是原函数的曲线，而绿色是通过fit得到的多项式函数，而蓝色点则是使用neville方法计算的0.6处的值。虽然比较两种拟合方式比较接近，但是由于我们的数据不太符合多项式的特性，因此拟合还是出现了较大的偏差，因此，再次说明，多项式拟合并不一定能够构造出客观反映数据规律的多项式。这种方式得到的在0.6处的拟合值结果为2.2488997156639994，两者几乎具有完全相同的结果。因此，可知在0.6处的拟合值，其结果为2.2488997156640003。，该值就是在对应值为x时的拟合值。

如果给定N个点，然后要求一个多项式通过这N个点，最简单直接的方式是列出线性方程求解，N个点可以确定N个未知量，则所求的拟合多项式，其最高次幂就是（N-1）。这个公式具有很强烈的轮换性，和线性方程组的克莱姆法则很相似（其实本质上就是一样的），所以实现起来也是非常简单。可以看到三个图形是几乎完全重叠的，特别是使用拉格朗日法的插值多项式b，与原多项式a完全相同。当然，上述函数并未对输入进行检测，例如x，y的尺寸应该是一致的等。本篇介绍一下多项式插值中，拉格朗日法的原理及其实现。

常见的插值（拟合）方法有拉格朗日法，Neville迭代法，牛顿差商法等，还有特定场合下的Hermite插值、三次样条插值等，将在后续逐步展开介绍。继续说多项式的数值及拟合，这次主要讨论关于多项式拟合的函数fit。还要指定其域（domain），该值通常应该设定为[-1,1]。和原来的多项式相同。

除了进行值替换，多项式之间可以进行代数的基本运算，如加减乘除等，甚至还可以进行微分和积分等常见运算。

可以看到，求解是在复数域内完成的，这个也是代数基本定理的：n次复系数多项式方程在复数域内有且只有n个根（m重根算m个），这个定理最先由德国数学家罗特，达朗贝尔、欧拉和拉格朗日都先后尝试证明但不完善，但是在高斯的博士论文中，给出了第一个严格证明（不愧是数学王子），据说目前已经有两百多种证明方法。从数学角度来说，a和b这两个多项式应该是一样的，但是从代码运行的结果可见，b在输出时，最高次幂虽然系数为0，但是仍旧显示出来了，而且b的度和a的度也不一样。

对于一个在区间[a,b]上连续的函数f(x)，如果f(a)*f(b)

牛顿法速度快效率高，但是有个很大的问题是必须要事先知道函数的导数，有时候计算导数是一件费时费力的事情，例如上例中是预选定义好了导函数。为了克服这个问题，可以逆向思维，导数的定义之一通过计算微商，然后取极限的方式，可以将。可以看到，只需要4次就可以计算出误差在1e-6的水平。测试显示通过6次迭代也可以实现误差在1e-6的水平。继续求根的数值算法。这次介绍牛顿法和割线法。

这部分讨论一下矩阵的求逆和线性方程组的求解。如果两个矩阵A，B的矩阵乘积为单位矩阵，则互为逆矩阵，且有AB*BA=I。事实上这种方式由于计算量太大，通常不大会使用，但是这个法则确实很优雅。另一种方式是使用逆矩阵。其中A为矩阵，x和b为列向量。

在numpy中，ndarray表示向量和矩阵，是N-dimensional array（N维数组）的缩写，由多个具有相同类型和尺寸的元素组成的多维容器。

本次讨论一下函数的定义和使用。

本篇介绍了均分向量的创建和使用，包括自定义实现和现成的软件包

介绍通过Python进行数值计算的基本运算和库