http://scz.617.cn:8/misc/202112152114.txt
这篇是为数学爱好者续写的。微软资深程序员Larry Osterman写过一篇blog
My son is SUCH a geek (in a good way) - [2007-11-07]
文中用"连续差分法"求解高阶等差数列的通项公式。我展开一下，做点数学科普。
先介绍下科普背景。这些年小学阶段有一些"班"，这些班有一些所谓的数学测试，给定一个数字序列，问下一个数字应该是什么。一般出题方认定只有一个"标准答案"，但从数学角度看，这种题有无穷多种逻辑自洽的答案，如果真是数学题，每给出一个答案就要加一些约束条件。这是一种扯淡的门萨式测试，Leo Davidson喷得精彩，我引用一下：
这就是为什么我不能忍受那些门萨式测试，这些自作聪明的家伙认为，你应该找出最简单、最显而易见的答案。但简单性和显而易见性完全是主观的，与数学几乎无关。
Larry Osterman的儿子Daniel在他五年级时遇到过这样一个数字序列
79 561 2279 6445 14679 29009
他当时的数学老师Bob Whittemore教了他一种算法，用于寻找该数字序列的"规律"，套用中学数学术语，即根据已知条件求解该数列通项公式。考虑到五年级这个前提，没有其他约束条件是可以理解的。
求解过程大致如下
x f(x)
1 79
2 561 482
3 2279 1718 1236
4 6445 4166 2448 1212
5 14679 8234 4068 1620 408
6 29009 14330 6096 2028 408
n=1 n=2 n=3 n=4
第1列是x，第2列是f(x)，认为它是一个任意数列。第3列是第2列相邻两项的差形成的新数列，依次类推。将第3列视为原数列(第2列)的第1次迭代求差，直至第n次迭代求差后出现等差数列，停止迭代求差，表示f(x)的最高次数是n，x^n的系数是"稳定差"除以n!。
本例n=4，稳定差等于408，408/4!=408/24=17，f(x)的最高次项是17*x^4。
现在将原数列(第2列)各项减去17*x^4，得到一个新数列，用同样的办法求解新数列的最高次项。显然，新数列经过3次迭代求差会出现稳定差。
x f(x)-17*x^4
1 79-17=62
2 561-17*16=289 227
3 2279-17*81=902 613 386
4 6445-17*256=2093 1191 578 192
5 14679-17*625=4054 1961 770 192
6 29009-17*1296=6977 2923 962 192
n=1 n=2 n=3
此次n=3，稳定差等于192，192/3!=192/6=32，"f(x)-17x^4"的最高次项是32x^3。
x f(x)-17*x^4-32*x^3
1 62-32=30
2 289-32*8=33 3
3 902-32*27=38 5 2
4 2093-32*64=45 7 2
5 4054-32*125=54 9 2
6 6977-32*216=65 11 2
n=1 n=2
此次n=2，稳定差等于2，2/2!=1，"f(x)-17x^4-32x^3"的最高次项是1*x^2。
x f(x)-17*x^4-32*x^3-x^2
1 30-1=29
2 33-4=29
3 38-9=29
4 45-16=29
5 54-25=29
6 65-36=29
现在有
f(x)-17*x^4-32*x^3-x^2=29
求得数列通项公式
f(x)=17*x^4+32*x^3+x^2+29
上述算法的数学术语是"连续差分法"，英文是"Method of Successive Differences"。
Daniel的轶事到此就结束了，后面是其他展开。
f(1)到f(6)都符合已知值
f(1)=171^4+321^3+1^2+29=79
f(6)=176^4+326^3+6^2+29=29009
由于有数列通项公式，可以接着求f(7)
f(7)=177^4+327^3+7^2+29=51871
相当于回答"下一个数字是什么"。但是，你再好好想想"连续差分法"的求解过程，完全可以给出一个与f(7)不同的数字，然后用"连续差分法"求解"新数列通项公式"，假设是g(x)，g(1)到g(6)与f(1)到f(6)相符，但g(7)不等于f(7)。g(7)这个答案逻辑自洽，从数学上讲，你不能说g(7)不是答案的一种。
随手在百度上搜了一个类似的题，问，已知数列
1 5 18 59 184
下一个数字应该是什么？再次演示"连续差分法"
x f(x)
1 1
2 5 4
3 18 13 9
4 59 41 28 19
5 184 125 84 56 37
n=1 n=2 n=3 n=4
第1列是x，第2列是f(x)，认为它是一个任意数列。第3列是第2列相邻两项的差形成的新数列，依次类推。将第3列视为原数列(第2列)的第1次迭代求差，直至第n次迭代求差后出现等差数列，停止迭代求差，表示f(x)的最高次数是n，x^n的系数是"稳定差"除以n!。
本例n=4，稳定差等于37，37/4!=37/24，f(x)的最高次项是37/24*x^4。
现在将原数列(第2列)各项减去37/24*x^4，得到一个新数列，用同样的办法求解新数列的最高次项。显然，新数列经过3次迭代求差会出现稳定差。
x f(x)-37/24*x^4
1 1-37/24=-13/24
2 5-37/24*16=-472/24
3 18-37/24*81=-2565/24
4 59-37/24*256=-8056/24
5 184-37/24*625=-18709/24
为将来计算方便，令
g(x)=-4!*3!*2!*1!*(f(x)-37/24*x^4)
g(x)=-24*6*2*(f(x)-37/24*x^4)
g(x)=-288*(f(x)-37/24*x^4)
为什么要乘以"4!3!2!1!"呢？因为"连续差分法"确定各次项系数时要除以n!，我不想约分来约分去的，乘以"4!3!2!1!"可以确保后续计算不再出现分数，等求出g(x)的多项式表达，但做简单变换求f(x)。
有
x g(x)
1 13/24*288=156
2 472/24*288=5664 5508
3 2565/24*288=30780 25116 19608
4 8056/24*288=96672 65892 40776 21168
5 18709/24*288=224508 127836 61944 21168
n=1 n=2 n=3
此次n=3，稳定差等于21168，21168/3!=21168/6=3528，g(x)的最高次项是3528*x^3。
x g(x)-3528*x^3
1 156-3528=-3372
2 5664-3528*8=-22560 -19188
3 30780-3528*27=-64476 -41916 -22728
4 96672-3528*64=-129120 -64644 -22728
5 224508-3528*125=-216492 -87372 -22728
n=1 n=2
此次n=2，稳定差等于-22728，-22728/2!=-22728/2=-11364，"g(x)-3528x^3"的最高次项是-11364x^2。
x g(x)-3528*x^3+11364*x^2
1 -3372+11364=7992
2 -22560+11364*4=22896 14904
3 -64476+11364*9=37800 14904
4 -129120+11364*16=52704 14904
5 -216492+11364*25=67608 14904
n=1
此次n=1，稳定差等于14904，14904/1!=14904，"g(x)-3528x^3+11364x^2"的最高次项是14904*x。
x g(x)-3528*x^3+11364*x^2-14904*x
1 7992-14904=-6912
2 22896-14904*2=-6912
3 37800-14904*3=-6912
4 52704-14904*4=-6912
5 67608-14904*5=-6912
现在有
g(x)-3528*x^3+11364*x^2-14904*x=-6912g(x)=3528*x^3-11364*x^2+14904*x-6912
g(x)=-288*(f(x)-37/24*x^4)
f(x)-37/24*x^4=-3528/288*x^3+11364/288*x^2-14904/288*x+6912/288
求得数列通项公式
f(x)=37/24*x^4-49/4*x^3+947/24*x^2-211/4*x+24
上式不方便计算，换一个
f(x)=(444*x^4-3528*x^3+11364*x^2-14904*x+6912)/288
验证
f(1)=(4441^4-35281^3+113641^2-149041+6912)/288=1
f(5)=(4445^4-35285^3+113645^2-149045+6912)/288=184
求f(6)
f(6)=(4446^4-35286^3+113646^2-149046+6912)/288=486
486只是答案之一，但不是唯一答案，有无穷多种答案。
接下来科普一下"高阶等差数列"，不做严格数学论述，拣要点讲讲。假设原数列是这种
79 561 2279 6445 14679 29009
或
1 5 18 59 184
"连续差分法"图示中"n=1"的列称为原数列的"1阶差数列"，"n=4"的列称为原数列的"4阶差数列"，假设"n=r"时首次出现"非零稳定差"，则称原数列为"r阶等差数列"。中学阶段所学的非常数等差数列就是"1阶等差数列"(不是"1阶差数列")，常数序列可以定义成"0阶等差数列"。
我国中学阶段常规教学会讲"1阶等差数列"通项公式，但不会讲如何求解高阶等差数列通项公式，事实上有成熟的定理，不需要显式动用"连续差分法"。
设{an}是r阶等差数列，其各阶等差数列首项为d1、d2、…、dr，则数列的通项公式是
an=a1+C(n-1,1)*d1+C(n-1,2)*d2+...+C(n-1,r)*dr
其中
C(n,m)=n!/m!/(n-m)!
即组合公式。只看这个公式可能有些懵逼，举个例子。
x f(x)
1 1
2 5 4
3 18 13 9
4 59 41 28 19
5 184 125 84 56 37
n=1 n=2 n=3 n=4n=4
a1=1
d1=4
d2=9
d3=19
d4=37
an=1+C(n-1,1)*4+C(n-1,2)*9+C(n-1,3)*19+C(n-1,4)*37
或者写成
f(x)=1+C(x-1,1)*4+C(x-1,2)*9+C(x-1,3)*19+C(x-1,4)*37
验证一下
f(5)=1+C(4,1)*4+C(4,2)*9+C(4,3)*19+C(4,4)*37
=1+4!/1!/3!*4+4!/2!/2!*9+4!/3!/1!*19+4!/0!/4!*37
=1+24/1/6*4+24/2/2*9+24/6/1*19+24/1/24*37
=184f(6)=1+C(5,1)*4+C(5,2)*9+C(5,3)*19+C(5,4)*37
=1+5!/1!/4!*4+5!/2!/3!*9+5!/3!/2!*19+5!/4!/1!*37
=1+120/1/24*4+120/2/6*9+120/6/2*19+120/24/1/*37
=486
r阶等差数列{an}的通项公式一定是关于正整数n的r次多项式，事实上
f(x)=1+C(x-1,1)*4+C(x-1,2)*9+C(x-1,3)*19+C(x-1,4)*37
=(444*x^4-3528*x^3+11364*x^2-14904*x+6912)/288
本文直接引用了数学结果，没有给出数学证明，有兴趣者去看相关代数书籍好了。
数值分析中通过已知离散数据求未知数据的过程或方法称为插值法。用解析几何的绘图来看，已知若干离散点坐标(xi,yi)，求一条曲线，同时经过所有已知离散点。显然，有无数条曲线经过所有已知离散点。这更直观地表明"下一个数字是什么"，有无穷多的答案。每条曲线都对应某一种插值法，最扯淡的当然是线性插值法，将所有已知离散点用折线连起来，这种插值法完全没有意义。
关于插值法，可以搜多项式插值法、牛顿插值法、拉格朗日插值法等关键字。我不是数学系的，数值分析这块白瞎，这些玩意儿的科普我就别丢人现眼了。但我会讲数学八卦啊，各位坐好了听。
尊重历史的话，牛顿插值法应该称为"格雷戈里-牛顿公式"，格雷戈里和牛顿各自独立给出该公式。但牛B顿这位，大家都认识的吧，神特么废话，属于"天不生牛顿、万古如长夜"的主儿，他在他那套想不知道也难的《自然哲学的数学原理》第三卷引理五里发表了牛顿插值法，后面的事大家都知道了。这个故事告诉我们，跟牛B顿抢C位，你得是莱布尼茨这种级数的，才好平分秋色。要么降一档，也得是胡克这种级数的，让牛B顿反讽地说句，如果我能看得更远一点的话，是因为我站在巨人的肩膀上。一般人真抢不过牛B顿。
格雷戈里(James Gregory)其实也很厉害。他小时候跟母亲学习几何学，妥妥的学霸家庭。应该是早期发明反射望远镜的人，后来反射望远镜以他的名字命名。这哥们是个狂热天文爱好者，由于长期使用望远镜进行天文观测而导致眼疲劳，结果眼瞎了。后来在爱丁堡去世。嗯，女魔头在那里读过书。
拉格朗日插值法最早被一个英国数学家发现，四年后由欧拉再次发现，又过了十二年，拉格朗日在他的《师范学校数学基础教程》中发表该插值法。最后出场的抢到了C位，该插值法被命名为拉格朗日插值法。话说能抢过欧拉，这主要凭运气，虽然拉格朗日实力强悍，但欧拉是那么好相与的？欧拉是谁？数学史上用手数得过来的主儿之一。
当然，我对拉格朗日很崇拜。拉格朗日平动点L4、L5的理论测算及科学运用，完美展示了自然哲学的数学原理，一种科学的奇迹。
过去是戏不够、烟来凑，我这是典型的没有数学天赋的货，只能东拉西扯些有的没的。
如有侵权请联系:admin#unsafe.sh