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的理论测算及科学运用,完美展示了自然哲学的数学原理,一种科学的奇迹。
过去是戏不够、烟来凑,我这是典型的没有数学天赋的货,只能东拉西扯些有的没的。