问题的起源

问题的起源是我去年在公司研究的一个问题。问题的背景是研究主权违约风险对外汇波动率的影响，在数学上可以化归为带跳的Brownian运动的分布问题。总之，在一系列复杂的运算之后（很明显，绝大多数的运算都不是我做的……），问题被转化为求解一个无穷重积分： 此处的有一个很好的性质： 当然在现实中并不能计算无穷重的积分，但好在积分的重数越多只是会让结果更精确。利用scipy包的三重积分，已经可以得到不错的结果——只可惜距离我们的理想值还存在差距。然而scipy包的重积分到三重积分就为止了，更高重积分就必须自己动手设计算法了。

单重积分的数值方法

我们自然是要牢记祖训「不要重复造轮子」，所以我们看看目前的工具箱里都有什么。首先我们看看单重积分有哪些数值方法。

从单重积分的几何意义上，我们可以很容易地得到两种朴素的近似方法：

  1. 中点法
  2. 梯形法

这两种方法看起来差不多，但是实际不然。我们可以证明如下定理：

定理：对于上的可导凸函数，中点法的误差不大于梯形法误差。

要证明这个定理，我们需要先证明一个引理。这个引理被称为Hermite-Hadamard（埃尔米特-阿达马）不等式（以下简称HH不等式）。

引理：对于上的可导凸函数，以下不等式成立：

证明： 令，其中，则 其中由Jensen不等式给出。

另一方面，由于是上的可导凸函数， 两边积分得

现在我们证明原命题： 令，由 代入HH不等式右边部分，可得 整理得

现在我们知道了对于凸函数，中点法比梯形法精确，但是我们还希望知道一个问题——到底有多精确？我们接下来就回答这个问题。

定理（中点法的精度）：设在上可求2次导，则使得： 证明： 由Taylor展开式： 两边积分得：

我们不加证明（证明可参考Mercer, P. R. (2014). More calculus of a single variable. Berlin: Springer.）地指出，梯形法的精度是：

☡： 此处我们并没有要求是上的凸函数。

我们注意到，虽然两个误差项的系数不同，但都包含了这一项。因此，我们给出一个（不严格的）定义：如果某种数值积分的方法误差项可以写做，我们就称这种方法有阶代数精度。在这里，中点法和梯形法都具有2阶代数精度。

我们还注意到，中点法的误差是，梯形法的误差是，那么能不能把两种方法「匀一匀」，创造出一种误差更小的方法呢？（虽然这种想法很危险，因为两边的是不一样的）

——还真的可以。我们通过组合的中点法和的梯形法得到了一种新的方法：

这种方法被称为Simpson法。（显然，这种方法和日后成为预言纪录片的『辛普森一家』没有任何关系……）

那么，Simpson法真的有更好吗？下面我们来尝试计算Simpson法的误差。首先，我们假设在上充分光滑（就是说能求很多阶导数，具体求多少阶我们还不知道），令并且构造一个新的辅助函数（别问我这是怎么来的，问就是『辛普森一家』预言出来的——开玩笑的，稍后揭晓） 不知道这东西可以干什么，于是我们求导试试看： 且 的形式看起来很不错，所以我们再做一些操作。

由Lagrange中值定理，，使得： 这样，我们就可以把写成

令，连续运用三次Cauchy中值定理可知，，使得： 故 而： （现在知道是怎么被构造出来了的吧）。 由此可知，Simpson法的代数精度是3阶。——通过「匀一匀」，还真的把精度提高了。

——还能更给力一些吗？其实可以。如果我们把区间均匀地分成n份，即令，，再对这个点进行多项式拟合（没错，就是传说中的那个Lagrange插值公式——那个可以让你在任何数列找规律的题目中都回答114514的神奇公式），最后再对这个多项式积分（多项式积分是由确切的公式的），理论上给出的点越多精度越高。我们可以证明有个点的方法就能达到阶代数精度。

现在我们通灵一下，直接给出5个点和7个点的情况（计算过程非常无聊，有兴趣的读者可自行验算w） （建议检查一下计算这些鬼公式的时候的贝拉的精神状态）

我们注意到，前面给出的阶代数精度误差公式都有，显然，当区间很长时这一项也会爆炸性增长，这是我们不想看到的。因此，我们可以将区间分成份，在每一个小区间上运用这些求积公式。令，，那么

代入中点公式： 代入Simpson公式： 通过这种方法，我们可以将阶代数精度误差公式中的缩小到，当区间分得够细，我们就可以得到更精确的结果。

重积分的数值方法

这里我们使用一个最朴素的想法：将重积分一层一层地用数值公式「剥开」。以中点公式为例： 再以Simpson公式为例： 这个展开式已经非常麻烦了，但是如果我们有性质 就可以简化为： 进一步地，如果有： 那么 其中 同样地，对于5个点和7个点的求积公式，我们也可以写出类似的公式：

我们也可以在这里运用将区间分成N份的方法，利用这种方法的中点公式： 同理： 代入Simpson公式：

最后的问题——时间复杂度

从上面我们可以看出，如果直接计算这个和式，它的时间复杂度是，也就是说如果取得很大的话，这个算法是很慢的。那么我们有没有办法在算法上改进呢？

我们注意到，只要有相同的和，那么最后的函数的值就是一样的，所以我们只需要考虑：有多少种可能的结果，使得，其中

这是一个典型的拆分数问题。

最初，在某位猫猫的提醒下，我发现将一个数拆分为个正整数等价于将个相同的球排列成一排，向其中的个空隙中插入个隔板，总共可能的情况。这就是我们高中学排列组合时典型的隔板法问题，结果是

后来经过验算，我发现这个结果并不对。原因是这并不是一个普通的拆分数问题，而还有一个附加条件：拆出来的最大数不能超过。

于是我们换一种思考方式：利用多项式。我们注意到，将拆分为个正整数，其中最大的数不能超过等价于求多项式中的系数，进一步地，等价于求多项式中的系数。

这里我们就需要用到一般母函数的理论了。注意到： 我们令，那么我们可以得到的系数是。 问题被完美解决。