植树节晚间，和和王JM谈到数值计算方法这节课。然后然后他说“我们有数值分析“。我问“数值分析是什么”，“数值计算方法plus”。之后为展示该门课为什么是Plus，他给我看了一个颇为有趣的题目：

Produce a table of the values of the series

${\sum }_{k=1}^{+\mathrm{\infty }}\frac{1}{k(k+x)}$
\sum\limits_{k=1}^{+\infty}\frac{1}{k(k+x)}

for the 3001 values of x, x= 0.000, 0.001, 0.002, …, 3.000. All entries of the table must have an absolute error less than 0.5e-10 (12 digits of precision). This problem is based on a problem from Hamming (1962), when mainframes were very slow by today’s microcomputer standards.

跟着王JM的思路，我以“Numerical Summation of A Series”为标题试着搜了一下这道题。果不其然，出现了一遛的CSDN……

其中一种非常典型的方法（采用裂项法）如下：

需要用公式将x不是整数的情况收敛到整数…

(不要问我为什么不用CSDN。我不知道我这个偏见是从哪里来的，但是我就是不用。)

至于为什么要这么做，该文档没有说。因为已经是深夜，直接问了王JM，他也说不出理论上的所以然，只是说如果硬计算，收敛速度不仅很慢，而且不准；而裂项又不好裂，通式并不唾手可得。不过，确实可以通过化成整数项的方式巧妙化简（真的简了吗？如简）。

但重点是，这个题目中提到了一个很有意思的人：

…This problem is based on a problem from Hamming(1962), when mainframes were…

这人是谁？

Richard Wesley Hamming (February 11, 1915 – January 7, 1998) was an American mathematician whose work had many implications for computer engineering and telecommunications. His contributions include the Hamming code (which makes use of a Hamming matrix), the Hamming window, Hamming numbers, sphere-packing (or Hamming bound), Hamming graph concepts, and the Hamming distance.

哥们，这么重量级的人物啊！！！汉明码创始人啊！！！作为计算机学院的人你怎么能忽略这个信息呢！！！！

令人寒心的是，搜遍全网关于该问题的解法，所有人无一例外，全部埋头苦算，硬算完了就没下文了，连解释都懒得解释，也没有好奇地去找找Hamming到底有没有写过这个玩意，甚至没有形成一片报告文档。这怎么行？你都不知道这道题来源于什么工程环境，你怎么就上去硬干了呢？这不工程师啊。

好吧，今天彻底成不眠之夜了。我非得挖出这道题来。

第一步，尝试搜索

本来如果Hamming在世，这个问题一封邮件就可以解决（且不论对方理不理我）。可惜生晚了。

我的第一想法是利用谷歌学术，穷举Hamming在1960年到1970年写过的所有文献。原因很简单，作者声称这个序列曾经被他在1962年开过光。但是我翻遍了Hamming在1960-1970写过的所有文献，没有。所有有关该序列的内容都一致指向2001年出现的ZOJ Problem Set，也就是该死的浙大自己的习题集中。

搜寻了两个小时后，我发现，最有可能的是，该例子作为某物的引证，出现在祖师爷的 Numerical Methods for scientists and engineers(1962) 中。这本书有足足七百多页，如果换白天，我可能会选择直接放弃。但是那是晚上，我非常不愿意把问题隔夜，不然会睡不着。

粗略二分搜寻了两遍，无果。实在撑不住，发了个朋友圈（我拒绝熬夜喝咖啡），发现张YH点了个赞。抱着哦，原来有人和我一样熬夜啊的这种心情，重整旗鼓，去目录上找。一个子条目引起了我的注意：

第十二章是无穷级数啊。看来这个地方大概会有。

功夫不负有心人，经过十分钟的仔细阅读，我发现了这个：

原来如此！这下我们就清楚这个无穷级数的庐山真面目了。这个东西其实就是伪装的非常好的Digamma函数(双伽玛函数。如果你不记得这玩意是什么的话，还记得Gamma函数指的是阶乘对实数域的拓展吗？双伽玛就是把他对数一下求个导)。其递归关系如下（不包括整数。整数的运算很简单，有另一套无穷级数，这里不讨论）：

$\psi (1+z)=-\gamma +{\sum }_{k=1}^{\mathrm{\infty }}(\frac{1}{k}-\frac{1}{z+k})$
\psi(1+z) = -\gamma + \sum\limits_{k = 1}^{\infty}(\frac{1}{k}-\frac{1}{z+k})

其中，γ是所谓的Euler–Mascheroni constant，欧拉-马斯刻若尼常数。具体细节可以在这里阅读。我们要的那玩意和他区别无非是：

${\sum }_{k=1}^{+\mathrm{\infty }}\frac{1}{k(k+x)}=\frac{1}{x}\left(\psi \right(x+1)+\gamma )$
\sum\limits_{k=1}^{+\infty}\frac{1}{k(k+x)} = \frac{1}{x}(\psi(x+1)+\gamma)

第二步，一种很简单，很简单的算法

大家如果熟悉Matlab，就应当会知道，这玩意是可以直接调用函数计算的。

和实际的真值完全没有误差。

Matlab的psi函数是什么原理？我打天落地查不到。但是通过Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables中的描述，我们可以猜测这东西要么是暴力打表的（事先存储[0,1]的一张表，用递推公式算。如果精度不够，那就再放大以后，再按照递推公式算！！），要么是通过任一种展开方式干的。原题用的是C语言做的，当然不能直接调用Matlab的库（这就是为什么计院的这门课不允许你用Matlab……控院的数值计算方法看起来还是挺文明的）。不过，我们可以站在巨人的肩膀上，看看别人是怎么做的。

我们来看一个几乎比我爸爸都老的论文：Bernardo, José M. (1976). “Algorithm AS 103 psi(digamma function) computation”。

1976年。按照当时的微处理器速度，累加出一大堆项是一个很烦的过程，更不用说打算加到一个很大的值让他收敛了。但是，谁说我们一定要用挫爆了的欧拉递推公式来算？

如果你的记忆力超群（比如说一些数院的天才——我显然不在列，我是控院的），那么你就会知道，有一种叫做Stirling’s approximation的东西，曾经被用在估计x值比较大的伽玛函数中。这东西可不可以用于估计双伽玛函数？显然是可以的。

$\psi (1+x)=ln\left(x\right)-\frac{1}{2x}-{\sum }_{j=1}^{\mathrm{\infty }}\frac{B_{2j}}{2j{x}^{2j}}$
\psi(1+x) = ln(x) - \frac{1}{2x} - \sum\limits_{j=1}^{\infty}\frac{B_{2j}}{2jx^{2j}}

这个级数对于任意x都收敛。其中，B代表的是伯努利数。

在实际应用中，我们只需要取这个级数前面的几项，比如，忽略掉7次以上的项，应该就足够了。但是这么干了之后，虽然说x值越大越精确，但对比较小的x而言（在这里是8.5左右以下），这个级数的符合程度就差到爆表了。 但是这怎么要紧？x小的时候的值可以用大的时候的值来推算啊。还记得这个吗：

$\psi (x+1)=\psi \left(x\right)+\frac{1}{x}$
\psi(x+1) = \psi(x) + \frac{1}{x}

那么，这堆小的不能再小的[0.001:3.000:0.001]，就可以这么来计算：

$\psi \left(x\right)=\psi (x+n)-{\sum }_{i=1}^{n-1}\frac{1}{x+i}$
\psi(x) = \psi(x+n) - \sum\limits_{ i = 1 }^{n-1}\frac{1}{x+i}

$\psi (x+n)\approx \ln (x+n)-\frac{1}{2(x+n)}-\frac{1}{12(x+n{)}^2}+\frac{1}{120(x+n{)}^4}-\frac{1}{252(x+n{)}^6}$
\psi(x+n) \approx \ln(x+n) - \frac{1}{2(x+n)} - \frac{1}{12(x+n)^2} + \frac{1}{120(x+n)^4} - \frac{1}{252(x+n)^6}

根据这个简单的递推，我们就可以来写程序了。我用的是Matlab，但是没有用到超过C语言范畴的东西。

1
2
3
4
5
6
7
8
9
10
11
12
13

function v = Digamma(y,n)
    x = y + n;
    sol = log(x)-1/(2*x)-1/(12*x^2)+1/(120*x^4)-1/(252*x^6);

    while(1)
        sol = sol - 1/( x - 1 );
        x = x - 1;
        if ( x > 1 && x < 2 ), break,end
    end
    v = sol;
end

结果展示如下。你可以自己去验证一下，符合地都很好。我没法都截取，这样太占据博客空间了。

诶，大就是好。那么n是越大越好吗？你此时已经嗅到了一丝不对的气息。我们试试看把n的值调大。

我去，有效位数怎么退化了？虽然说还是符合题目要求的十位，但是明显感觉到后面的玩意不对劲了。这样下去，之后一通计算猛如虎，误差变成二百五。

造成这样现象的原因很简单：我们采取这种方法的根本原因，就是要避开 带有和大整数一起玩的 很小的数的 项！但如果你把n的值设的太大，n和需要的x一对比，比如100000对比0.0001，这十三点东西又出现了！因此，n千万不能选的太大。经过不断测试，一百也许也太大了。区间在30到50是比较合理的，适用范围也比较广。

ok，最后为了满足题目条件，随便做个计算：

和题目答案一对比，丝毫不差。精度、速度还远远高于题目要求。何乐不为？

结语

双伽玛函数不止这种简单的计算方法，光维基百科上就能有近十种。我还有一大堆算法并没有拿出来讨论，因为他们的效果基本上和这种简单算法差不多，大同小异。不过也有很多算法很标新立异，比如这篇牛津大学的文献,用了一些令人瞠目结舌的方法进行计算。我觉得这彻底超纲了，就没有讨论，有兴趣，有能力的时候再自己做吧。

这个故事告诉我们，有问题别特么就想着硬算，先看看这个东西的本质是什么！

ps 关于本篇采用的算法：这篇1976年的论文居然是用Fortran写的。我不会这种语言，他写的也实别是烦逼醪糟的。最后自己写的玩意可能和他搞的有出入，烦请见谅。

reference

1.Hamming, R. (1962). Numerical methods for scientists and engineers. Dover Publications.
2.Abramowitz, M., & Stegun, I. A. (Eds.). (1968). Handbook of mathematical functions with formulas, graphs, and mathematical tables (Vol. 55). US Government printing office.
3.Bernardo, José M. (1976). Algorithm AS 103 psi(digamma function) computation.
4.Beal, Matthew J. (2003). Variational Algorithms for Approximate Bayesian Inference.

救命，忙不过来了

Assignments for the Numerical Methods of Zhejiang University by Dr.Red. Update regularly every Sunday.

浙大数值计算方法作业。俺写的。不要问我为什么用英文，好看。

Computer Assignment #1

普通的matlab使用练习。很挫，当时对latex的该模版还不是很熟悉，不建议观看。

Computer Assignment #2

虽然是小作业，还是用很长的时间改进了排版。应该花了一天左右，同时学习了Matlab如何画图。这里的精华在于第三种方法实际上是用牛顿法的特例（几千年前的海伦法）搞出来的，严格来说已经不是原先的那种方法了（二次收敛，而非线性收敛。这就是为什么迭代次数骤降了一倍），起码我拒绝承认这还是原来的那种纯粹的不动点迭代法。这个东西需要一定的观察才能看出来，不知道出出来是不是坑人的。

Computer Assignment #3

非常普通的矩阵求解。注意这个矩阵因为行列式接近于0，是病态的。另外对代码排版的探索也是在这一版中逐渐明确的

Computer Assignment #4

插值。可以感受到非常明显的Romberg效应，因此这个例子其实并不怎么适合lagrange插值，反而倒是真·线性插值适合一些。

Assignment #1

复习了一遍积分的递推公式。这里很有意思的点在于，递推需要考虑算法的收敛性。一开始正着递推误差会阶乘式爆炸性变大，而倒着推就会减小，所以采用倒推法。不过这需要超前的Newton-Cotes法等知识，不建议使用。

Assignment #2

求根的非常多的方法。很简单，没什么好说的。值得注意的是不动点法如果不经过转换，会直接飞出去。所以千万不要直接上。

Assignment #3

这篇年代有些久远了，技术细节基本上忘光了。要注意的大概是这道题的中文翻译表述有很严重的问题。原题写的“待萃取物的组分Yin”实际上指的是萃余相中溶质A的浓度，“萃取剂的组分Xin”是萃取相中溶质A的浓度，这两个描述的是同一个物质在不同萃取相中的组分差别，不是萃取剂自己的差别。不然可能会觉得“萃取剂”在“被萃取的过程中”正在“越来越少”，这显然错误。和赵豫红老师沟通过，她大概会改一下这道题的表述。

Assignment #4

非常美的一篇小报告。Gauss-Newton法在非线性回归中无往不利，百战百胜，唯一的问题是不能用在最简单的线性回归中。至于分析，嗯……马后炮而已，不要管。

Assignment #5

Gauss-Legendre方法很惊艳，但是确定Gauss点的方法更加惊艳。可以在维基百科上搜索GaussQuadrature的确定方法（包括Legendre多项式），让人大开眼界。

P.S. 很有趣的地方在于，当时问的Chatgpt，Chatgpt给出了一个让人摸不着头脑的副对角数值计算方法。后来找到了很冷门的文献，发现这正好是正确结果的计算式子的倒数表达形式……令人叹服AI的信息知会能力。

End-term Assignment

大作业做的计算化学，浅要地用休克尔规则解决了一下之前在做化学竞赛的时候困扰我的各种问题。反正有这个计算方法以后，对于大的有机碳骨架共轭体系应该不是什么大问题的，无非还要解决各种杂原子的微扰动。不过那就是别的问题了，作为搞控制的我也不用太拘泥于此。了却心愿而已。

Numerical Methods Using MATLAB by John H.Mathews.

你好！Red残破的小站点从今天开始就建立起来了。现在这个该死东西还没有经过整理和打扫，不过之后应该会有空就修整一下的。

这个站点会拿来干什么？我不知道。现在我还在学习html、Markdown之类的内容。现在大概只是用来存放一些我觉得很有趣的东西，防止换电脑以后丢掉的吧。目前而言还没有买域名的计划，github自带的免费域名应该就够了。

据说伟大的敲电脑的都会建立自己的个人博客。我懒，我先这样了。有缘再写吧。

/……施工中……/