缓慢增长数列的逼近形式

问题

数学星空于2015年底提问 :
对于a(n+1)=a(n)+1a(n)2,a(1)=1a(n+1)=a(n)+\frac{1}{a(n)^2},a(1)=1, 我们可以利用微分方程求解a(n)a(n)的第一阶项. 设a(n)=ya(n)=ya(n+1)a(n)=a(n)a(n+1)-a(n)=a^{\prime}(n) 则有yy2=1y^{\prime}y^2=1 求解得到:a(n)=y(3n)13a(n)=y \approx (3n)^{\frac{1}{3}} 接下来,我们可设 (a(n))3=3n+a+bln(n)+cn+dn2+en3+fn4+(a(n))^3=3n+a+b\ln(n)+\frac{c}{n}+\frac{d}{n^2}+\frac{e}{n^3}+\frac{f}{n^4}+\dots 本主题的任务是如何求解a,b,c,d,ea,b,c,d,e?

或者利用数学计算软件或者数学理论给出更精确的表达式?

注:Kuing(郭子伟)已证明(P111,kuingluing20151222)

(3n)13<a(n)<(3n)13+(3n)13(3n)^{\frac 13}<a(n)<(3n)^{\frac13}+(3n)^{-\frac13}

并指出(彩色の夢∩o∩)给出下面结果,但并未给出证明过程:

(a(n))3=3n+2+ln(n+2)19n+6ln(5)+115(a(n))^3=3n+2+\ln(n+2)-\frac{1}{9n+6}-\ln(5)+\frac{1}{15}

最终我们通过努力,得出这个数列的无穷项展开形式。而且这种方法可以使用在很多类似的数列上,比如:
2010年初数学星空提问的一个数列不等式 也可以用类似方法分析:
若数列an{a_n}满足 :
a1=12,an+1=12(an2+1)a_1=\frac12 , a_{n+1}=\frac12(a_{n}^2+1)
则: 12n+2n2ln(n3)+417128n2an12n+5lnn+32n21-\frac2n+\frac2{n^2}\ln(\frac n3)+\frac{417}{128n^2}\le a_n\le 1-\frac 2n+\frac{5\ln n+3}{2n^2}

初步逼近

mathe很快给出如下分析方案 :
a(n+1)=a(n)+1a(n)2a(n+1)=a(n)+\frac{1}{a(n)^2}
得出
a(n+1)3=a(n)3+3+3a(n)3+1a(n)6a(n+1)^3=a(n)^3+3+\frac{3}{a(n)^3}+\frac{1}{a(n)^6}
归纳容易得出对于n2n\ge 2a(n)33n+2a(n)^3\ge 3n+2于是重新代入上面得出
a(n)38+3(n2)+3k=2n113k+2+k=2n11(3k+2)2a(n)^3\le 8+3(n-2)+3\sum_{k=2}^{n-1}\frac{1}{3k+2}+\sum_{k=2}^{n-1}\frac{1}{(3k+2)^2}
然后我们可以利用积分估算得出
a(n)38+3(n2)+ln(3(n1)+2)ln(8)+1813(n1)+2=3n+ln(3n1)13n1+18+2ln(8)<3n+2+ln(n)a(n)^3\le 8+3(n-2)+\ln(3(n-1)+2)-\ln(8)+\frac{1}{8}-\frac{1}{3(n-1)+2}=3n+\ln(3n-1)-\frac{1}{3n-1}+\frac{1}{8}+2-\ln(8)\lt 3n+2+\ln(n)
而在得到这个上界后,我们有可以利用它来估计下界得出
a(n)38+3(n2)+3k=2n113n+2+ln(n)a(n)^3 \ge 8+3(n-2) +3\sum_{k=2}^{n-1} \frac{1}{3n+2+\ln(n)}
而且容易看出后面求和式3k=2n113n+23\sum_{k=2}^{n-1}\frac{1}{3n+2}的差在nn趋向无穷收敛到一个常数
所以可以知道存在常数CC是的a(n)33n+ln(n)+Ca(n)^3 \ge 3n+\ln(n)+C
也就是我们可以非常轻松证明a(n)33nln(n)a(n)^3-3n-\ln(n)有界,进一步分析还可以得出其极限存在,但是求这个极限比较有难度
而数值计算Δ(n)=a(n)33nln(n)\Delta(n)=a(n)^3-3n-\ln(n),有
Δ(10)=1.220832\Delta(10) =1.220832
Δ(100)=1.151530\Delta(100)=1.151530
Δ(1000)=1.137657\Delta(1000)=1.137657
Δ(10000)=1.135573\Delta(10000)= 1.135573
Δ(100000)=1.135295\Delta(100000)=1.135295
Δ(1000000)=1.135261\Delta(1000000)= 1.135261
Δ(10000000)=1.135256\Delta(10000000)=1.135256
所以可以看出1#彩色の夢∩o∩)估计公式中常数项偏离比较大,倒是本楼开头估计是中对应常数项ln(3)+18+2ln(8)\ln(3)+\frac{1}{8}+2-\ln(8)只略微大了一些
同样要数值计算这个极限值,也可以先计算前面若干项,然后对于余下部分求和式和近似公式求差值,而差值部分误差可以非常容易估计,这样就可以得出极限值一个较好的估计范围了。

无穷展开

次日Buffalo给出了一个神奇的含ln(n)\ln(n)无穷展开式形式 :
首先令bn=an3b_n=a_n^3,得到bn+1=bn+3+3bn+1bn2b_{n+1}=b_n+3+\frac{3}{b_n}+\frac{1}{b_n^2}。渐近式前三项当然是3n+lnn+a3n+\ln n+a,但是后面的项可不是1n\frac{1}{n}的简单级数,而应该是i=1Pi(lnn,a)ni\sum_{i=1}^{\infty}\frac{P_i(\ln n,a)}{n^i},这里的Pi(lnn,a)P_i(\ln n,a)lnn\ln nii次多项式(事实上同时也是aaii次多项式)。用待定系数法,令fm(n)=3n+lnn+a+i=1mPi(lnn,a)nif_m(n)=3n+\ln n+a+\sum_{i=1}^{m}\frac{P_i(\ln n,a)}{n^i}让mathematica在n=n=\infty处展开fm(n+1)fm(n)33fm(n)1fm(n)2f_{m}(n+1)-f_m(n)-3-\frac{3}{f_m(n)}-\frac{1}{f_m(n)^2}O(1nm+1)O(\frac{1}{n^{m+1}})逐次得到Pi(lnn,a)P_i(\ln n,a)。计算发现Pi(lnn,a)P_i(\ln n,a)实际上很神奇地仅仅是是lnn+a\ln n+aii次多项式,而且系数有很明显的规律,可以把一部分项合并降低次数,目前做到bn3n+lnn+a+ln(1+lnn+a3n)56(3n+lnn+a)+ln(1+lnn+a3n)3n+lnn+a+i=2Qi(lnn+a)ni+2b_n\approx 3n+\ln n+a+\ln(1+\frac{\ln n+a}{3n})-\frac{5}{6(3n+\ln n+a)}+\frac{\ln(1+\frac{\ln n+a}{3n})}{3n+\ln n+a}+\sum_{i=2}^{\infty}\frac{Q_i(ln n+a)}{n^{i+2}},想进一步降低多项式次数遇到数列{270, 870, 1725, 2780, 4000, 5361},找不到规律。

前的结果用于快速计算任意初值的数列极限值已经足够了:已经知道bn3n+lnn+a+ln(1+lnn+a3n)56(3n+lnn+a)+ln(1+lnn+a3n)3n+lnn+a16n2+i=1Qi(lnn+a)ni+2b_n\approx 3n+\ln n+a+\ln(1+\frac{\ln n+a}{3n})−\frac{5}{6(3n+\ln n+a)}+\frac{\ln(1+\frac{\ln n+a}{3n})}{3n+\ln n+a}-\frac{1}{6n^2}+\sum_{i=1}^{\infty}\frac{Q_i(\ln n+a)}{n^{i+2}},从初值出发计算NN项,定义Δ(n)=bn3nlnn\Delta(n)=b_n-3n-\ln nΔ(N)\Delta(N)就是极限a的粗略值d0d_0,然后令di+1=Δ(N)ln(1+lnN+di3N)+56(3N+lnN+di)ln(1+lnN+di3N)3N+lnN+di+16N2d_{i+1}=\Delta(N)-\ln(1+\frac{\ln N+d_i}{3N})+\frac{5}{6(3N+\ln N+d_i)}-\frac{\ln(1+\frac{\ln N+d_i}{3N})}{3N+\ln N+d_i}+\frac{1}{6N^2},这样迭代数次直至结果不变,这个值和真实极限值的误差是O(lnNN3)O(\frac{\ln N}{N^3})。 费点功夫把展开式多计算几项就可以得到更好的逼近效果O((lnN)iNi+2)O(\frac{(\ln N)^i}{N^{i+2}})。 作为例子,选初值为1,计算100项,迭代几次后得到极限值为1.1352567349497837,计算10000项,迭代几次后得到极限值为1.1352558474127066,两者相差8.9×1078.9\times 10^{-7},小于估计的O(lnNN3)5×106O(\frac{\ln N}{N^3})\approx 5\times 10^{-6}

假设已经求得初值为b1b_1的数列的通式bn=f(n,b1)3n+lnn+hb1+R(n,a(b1))b_n=f(n,b_1)\approx 3n+\ln n+h_{b_1}+R(n,a(b_1)),现在来看对极限值hb1h_{b_1}我们能了解到什么程度。 递推方程是平移不变的,因此f(n,bi)=f(n,f(i,b1))=f(n+i1,b1)f(n,b_i)=f(n,f(i,b_1))=f(n+i-1,b_1),所以有hbi=limn(f(n,bi)3nlnn)=limnf(n+i1,b1)3nlnn=hb1+3i3h_{b_i}=\lim_{n\to\infty}(f(n,b_i)-3n-\ln n)=\lim_{n\to\infty}f(n+i-1,b_1)-3n-\ln n=h_{b_1}+3i-3,这是个严格的等式。由bi=f(i,b1)3i+lni+hb1+R(i,hb1)b_i=f(i,b_1)\approx 3i+\ln i+h_{b_1}+R(i,h_{b_1})可以迭代求得3ibilnbi3hb1+...3i\approx b_i-\ln\frac{b_i}{3}-h_{b_1}+...,最后得到hxxlnx33+56x+23x2+77108x3+133240x426695400x51676567x6+O(1x7)h_{x}\approx x-\ln \frac{x}{3}-3+\frac{5}{6x}+\frac{2}{3x^2}+\frac{77}{108x^3}+\frac{133}{240x^4}-\frac{2669}{5400x^5}-\frac{1676}{567x^6}+O(\frac{1}{x^7})。非常神奇的是尽管f(n,a)f(n,a)里含有巨量的lnn+a\ln n+a,可是后面的项里面不仅hb1h_{b_1}合情合理地消失了,连lnx\ln x也无影无踪。 对于很大的初值b1b_1,可以直接用这个表达式作为近似值。如果b1b_1不太大,可以先迭代i1i-1次得到较大的bib_i,算出hbih_{b_i},然后用等式hb1=hbi3i+3h_{b_1}=h_{b_i}-3i+3计算。

mathe对Buffalo的结果做进一步推导 :
对于bn+1=bn+3+3bn+1bn2b_{n+1}=b_n+3+\frac{3}{b_n}+\frac{1}{b_n^2},假设 bn=3n+ln(n)+a+k=1+Pk(ln(n)+a)nkb_n=3n+\ln(n)+a+\sum_{k=1}^{+\infty}\frac{P_k(\ln(n)+a)}{n^k}

我们先查看Pk(ln(n+1)+a)(n+1)k=Pk(ln(n)+a+ln(1+1n))nk×(1+1n)k\frac{P_k(\ln(n+1)+a)}{(n+1)^k}=\frac{P_k(\ln(n)+a+\ln(1+\frac{1}{n}))}{n^k\times (1+\frac{1}{n})^k}, 这个式子展开后是一个无穷级数,可以写成形如h=k+Qk,h(ln(n)+a)nh\sum_{h=k}^{+\infty}\frac{Q_{k,h}(\ln(n)+a)}{n^h}的级数,其中 每个Qk,h(.)Q_{k,h}(.)都是次数不超过k的多项式,而且同aa无关。
所以从形式上,这种形式的级数应该是可以用于这种递推式的,但是得出的结果不一定是收敛的,如同Γ\Gamma函数的Stirling级数就是不收敛的,而只是一种渐近式。
为了计算方便,我们可以将bn=3n+ln(n)+a+k=1+Pk(ln(n)+a)nkb_n=3n+\ln(n)+a+\sum_{k=1}^{+\infty}\frac{P_k(\ln(n)+a)}{n^k}中n用1y\frac{1}{y}代替,而ln(n)+a\ln(n)+aXX代替,写成形式级数
bn=3y+X+k=1+Pk(X)ykb_n=\frac{3}{y}+X+\sum_{k=1}^{+\infty}P_k(X)y^k
而其中如果我们将yy替换为y1+y\frac{y}{1+y},然后再将XX替换为X+ln(1+y)X+\ln(1+y),就可以得出bn+1=3(1+y)y+X+ln(1+y)+k=1+Pk(X+ln(1+y))yk(1+y)kb_{n+1}=\frac{3(1+y)}{y}+X+\ln(1+y)+\sum_{k=1}^{+\infty}\frac{P_k(X+\ln(1+y))y^k}{(1+y)^k} 然后代入原始递推式即可。
比如在Pari/Gp,我们可以输入

(16:47) gp &gt; 3/y+X+(b00+b01\times X)\times y+(b10+b11\times X+b12\times X^2)\times y^2+O(y^3)
%1 = 3\times y^-1 + X + (b01\times X + b00)\times y + (b12\times X^2 + b11\times X + b10)\times y^2 + O(y^3) //这个即$b_n$二阶近似
(16:49) gp &gt; subst(%1,y, y/(1+y))
%3 = 3\times y^-1 + (X + 3) + (b01\times X + b00)\times y + (b12\times X^2 + (-b01 + b11)\times X + (-b00 + b10))\times y^2 + O(y^3)
(16:51) gp &gt; subst(%3,X,X+log(1+y))
%4 = 3\times y^-1 + (X + 3) + (b01\times X + (b00 + 1))\times y + (b12\times X^2 + (-b01 + b11)\times X + (-b00 + (b01 + (b10 - 1/2))))\times y^2 + O(y^3)//这个为$b_{n+1}$的二阶近似
(16:52) gp &gt; %4-%1-3-3/%1-1/%1^2
%6 = ((-b01 + 1/3)\times X + (-b00 + (b01 - 11/18)))\times y^2 + O(y^3)

由此得出b01=1/3,b00=-5/18.
类似可以得出
bn=3n+ln(n)+a+13(ln(n)+a)518n+118(ln(n)+a)2+1154(ln(n)+a)16n2+181(ln(n)+a)3781(ln(n)+a)2+29162(ln(n)+a)1571458n3+1324(ln(n)+a)4+8243(ln(n)+a)3115972(ln(n)+a)2+122729(ln(n)+a)13327174960n4+...b_n=3n+\ln(n)+a+\frac{\frac{1}{3}(\ln(n)+a)-\frac{5}{18}}{n}+\frac{-\frac{1}{18}(\ln(n)+a)^2+\frac{11}{54}(\ln(n)+a)-\frac{1}{6}}{n^2}+\frac{\frac{1}{81}(\ln(n)+a)^3-\frac{7}{81}(\ln(n)+a)^2+\frac{29}{162}(\ln(n)+a)-\frac{157}{1458}}{n^3}+\frac{-\frac{1}{324}(\ln(n)+a)^4+\frac{8}{243}(\ln(n)+a)^3-\frac{115}{972}(\ln(n)+a)^2+\frac{122}{729}(\ln(n)+a)-\frac{13327}{174960}}{n^4}+...

公式:
k=1+Pk(X+ln(1+y))yk(1+y)kk=1+Pk(X)yk=3bn+1bn2\sum_{k=1}^{+\infty}\frac{P_k(X+\ln(1+y))y^k}{(1+y)^k}-\sum_{k=1}^{+\infty}P_k(X)y^k=\frac{3}{b_n}+\frac{1}{b_n^2} (3y+X+k=1+Pk(X)yk)2(k=1+Pk(X+h=1(1)h1yhh)yk(1+y)kk=1+Pk(X)yk)=3(3y+X+k=1+Pk(X)yk)+1(\frac{3}{y}+X+\sum_{k=1}^{+\infty}P_k(X)y^k)^2(\sum_{k=1}^{+\infty}\frac{P_k(X+\sum_{h=1}^{\infty}\frac{(-1)^{h-1}y^h}{h})y^k}{(1+y)^k}-\sum_{k=1}^{+\infty}P_k(X)y^k)=3(\frac{3}{y}+X+\sum_{k=1}^{+\infty}P_k(X)y^k)+1

上面表达式中常数a的值的计算是一个大问题,Buffalo最早给出了7位左右的精度,后来wayne 和mathe都尝试予以改进, Buffalo 随后给出了 30位精度的值1.135255847315503714195394347748。

快速计算方案

Buffalo随后表示:
绕开数列通式的渐近式直接快速求极限值的渐近展开式方法。
根据前面的分析,我们知道h(bi)=h(b1)+3i3h(b_i)=h(b_1)+3i-3,因为b2=b1+3+3b1+1b12b_2=b_1+3+\frac{3}{b_1}+\frac{1}{b_1^2},所以h(x)h(x)应该满足函数方程h(x+3+3x+1x2)=h(x)+3h(x+3+\frac{3}{x}+\frac{1}{x^2})=h(x)+3,从这个函数方程即可快速求得极限值a(x)a(x)的渐近式:前面的几项为
h(x)xlnx33+56x+23x2+77108x3+133240x426695400x51676567x6788279317520x7+7762807362880x8+88607258312247200x921193112111375x105126955274928226880x11+1317981395372314010796800x12+7041436522468181127498250880x13+973496293904381218979125x14+h(x)\approx x-\ln\frac{x}{3}-3+\frac{5}{6x}+\frac{2}{3x^2}+\frac{77}{108x^3}+\frac{133}{240x^4}-\frac{2669}{5400x^5}-\frac{1676}{567x^6}-\frac{788279}{317520x^7}+\frac{7762807}{362880x^8}+\frac{886072583}{12247200x^9}-\frac{21193112}{111375x^10}-\frac{51269552749}{28226880x^11}+\frac{13179813953723}{14010796800x^12}+\frac{7041436522468181}{127498250880x^13}+\frac{97349629390438}{1218979125x^14}+\cdots 除了系数可以乘以n(n+1)!3nn(n+1)!3^n化为整数实在是看不出规律,连符号都无法预测。

公式证明

mathe尝试对上面的计算过程给出比较严格的证明过程 :
f(x)=x+3+3x+1x2,f1(x)=f(x),fk+1(x)=f(fk(x))f(x)=x+3+\frac{3}{x}+\frac{1}{x^2},f_1(x)=f(x),f_{k+1}(x)=f(f_k(x))。那么fk(x)f_k(x)相当于保持数列b(n)b(n)递推关系不变,但是第零项b(0)b(0)改变成xx,那么第kk项就变为fk(x)f_k(x)了。
y=g(a,m)y=g(a,m)为使得limnfn(y)3(n+m)ln(n+m)=a\displaystyle\lim_{n\to\infty}f_n(y)-3(n+m)-\ln(n+m)=a的数而且y2y\ge 2
根据前面的分析过程,不管数列fn(y)f_n(y)中第零项yy是多少,极限limnfn(y)3nln(n)\displaystyle\lim_{n\to\infty}f_n(y)-3n-\ln(n)都是存在的,
于是对于给定的m极限limnfn(y)3(n+m)ln(n+m)\displaystyle\lim_{n\to\infty}f_n(y)-3(n+m)-\ln(n+m)也存在只是和前一个极限差了常数3m3m
于是y=g(a,m)y=g(a,m)表示如果保持数列b(n)b(n)的递推关系,其中第mmb(m)b(m)yy,那么limnb(n)3nln(n)\displaystyle\lim_{n\to\infty}b(n)-3n-\ln(n)的极限正好为aa
而比较有意思的是这时候的mm我们可以推广到任意一个实数而不仅仅是整数,从而将数列ba(m)b_a(m)推广为一个定义在正实数范围以mm为自变量的函数
于是容易得出g(a,m)=g(a+3d,md)g(a,m)=g(a+3d,m-d),而且对于我们以前定义的ba(m)=g(a,m)b_a(m)=g(a,m)
而根据前面分析应该有g(a,n)3n+ln(n)+a+k=1Pk(ln(n)+a)nkg(a,n)\approx 3n+\ln(n)+a+\sum_{k=1}^{\infty}\frac{P_k(\ln(n)+a)}{n^k}
由于g(a,n)=g(a3d,n+d)g(a,n)=g(a-3d,n+d),如果我们选择适当d使得ln(n+d)+a3d=0\ln(n+d)+a-3d=0,记其中m=n+dm=n+d,得到
那么结果可以简化为g(a,n)=g(a3d,m)3m+k=1Pk(0)mkg(a,n)=g(a-3d,m) \approx 3m+\sum_{k=1}^{\infty}\frac{P_k(0)}{m^k}.
由此,对于任意n,ba(n)=g(a,n)n,b_a(n)=g(a,n), 如果我们可以先利用上面公式反解出对应的m,就可以求出对应的d=mnd=m-n,然后利用a=3dlog(m)a=3d-log(m)得出aa.
于是问题就变成求渐进式s(x)3x+k=1Pk(0)xks(x)\approx 3x+\sum_{k=1}^{\infty}\frac{P_k(0)}{x^k},以及其逆函数.
然后我们假设当将n替换成n+1时,对应的mm变化为mm^{\prime},于是根据m的定义必然有mln(m)3=m+1ln(m)3m^{\prime}-\frac{\ln(m^{\prime})}{3}=m+1-\frac{\ln(m)}{3}
这里mm^{\prime}可以看成是mm的隐函数,显然可以看出在mm \to \infty时,m=m+1+o(1)m^{\prime} = m+1+o(1),代入递推式m=m+1+ln(mm)3m^{\prime}=m+1+\frac{\ln(\frac{m^{\prime}}{m})}{3},可以看出对于一般情况的f,这里m,mm^{\prime},m之间关系式直同f中常数项相关
得出m=m+1+13m+o(1m)m^{\prime} = m+1+\frac{1}{3m}+o(\frac{1}{m}),反复迭代可以得出m=J(m)=m+1+13m118m2154m3+7324m4+m^{\prime}=J(m)=m+1+\frac{1}{3m}-\frac{1}{18m^2}-\frac{1}{54m^3}+\frac{7}{324m^4}+\dots
系数依次为13,118,154,7324,314860,10729160,2833612360,621459270,\frac13,-\frac1{18},-\frac1{54},\frac7{324},-\frac{31}{4860},-\frac{107}{29160},\frac{2833}{612360},-\frac{621}{459270},\dots 然后再次利用公式ba(n+1)=ba(n)+3+3ba(n)+1ba(n)2b_a(n+1)=b_a(n)+3+\frac3{b_a(n)}+\frac1{b_a(n)^2},并且ba(n+1)=s(m),ba(n)=s(m)b_a(n+1)=s(m^{\prime}),b_a(n)=s(m),将m=J(m)m^{\prime}=J(m)代入,
也就是求s(J(m))=s(m)+3+3s(m)+1s(m)2s(J(m))=s(m)+3+\frac3{s(m)}+\frac1{s(m)^2}
比较系数,就可以得出Pk(0){P_k(0)}的值,
分别为518,16,1571458,13327174960,4441917873200,14533499330674400,77693360920832487200-\frac5{18},-\frac16,-\frac{157}{1458}, -\frac{13327}{174960},-\frac{444191}{7873200},-\frac{14533499}{330674400},-\frac{776933609}{20832487200}
从前面的结果来看都是负数而且绝对值小于1,也就是级数很可能是收敛的
于是H(x)=1s(1x)H(x)=\frac{1}{s(\frac1x)}满足H(0)=0而且展开式非常容易算出,于是H的逆函数在0的展开式也可以计算出来,由此可以得出s的逆函数的展开式,也是罗兰级数形式.
可以用来计算更高精度结构,然后结果再倒数一下可以得到s的逆函数形式,也就是可以用ba(n)b_a(n)直接计算对应的mm,最后可以推导出对应的aa
最后m的公式应该是13x+518x+12x2+239324x3+52276480x4+\frac13x+\frac5{18x}+\frac1{2x^2}+\frac{239}{324x^3}+\frac{5227}{6480x^4}+\dots.
前面得出m(x)的形式是m(x)=x3+k=1ukxkm(x)=\frac x3+\sum_{k=1}^{\infty}\frac{u_k}{x^k}。于是当x取bnb_n得出结果为m,而当x取bn+1b_{n+1}时得到结果为mm^{\prime},那么mm^{\prime}和m恰好有关系m=J(m)m^{\prime}=J(m).
也就是J(m(x))=m(x+3+3x+1x2)J(m(x))=m(x+3+\frac3x+\frac1{x^2}).
由此直接可以解出m(x)m(x)的系数.

在得到m(x)后可以求出对应的m再求出a.在得到a和给定项数n,那么可以先算m使得mln(m)3=n+a3m-\frac{\ln(m)}3=n+\frac a3,然后就可以用mxm\to x的函数s求出这一项的值.

现在比如选择b30=94.577173180254437058800379657718614759b_{30}= 94.577173180254437058800379657718614759
那么根据前4项可得m=31.5287182224435829596913315304969568m=31.5287182224435829596913315304969568,a=1.13525584723432304669623826092662724a=1.13525584723432304669623826092662724,精度在8位
如果b1000=3008.045412232225615137583957841337992049977740047149556926283b_{1000}=3008.045412232225615137583957841337992049977740047149556926283
那么根据m(x)前4项可得m=1002.681896477635973673360401688986574721554051620745012411842,a=1.1352558473155037114596749643830028099868785719581459235494458m=1002.681896477635973673360401688986574721554051620745012411842,a=1.1352558473155037114596749643830028099868785719581459235494458精度在17位,基本接近1x5\frac1{x^5}左右的精度.

而如果我们计算m(x)更多项数,就可以有更好的收敛精度,甚至对于较小的x也可以计算出高精度的a.

推广

wayne随后提出对问题进行一般化推广
此题目抽象出来,本质上是一个非线性函数在无穷迭代后的函数性态分析。 设x1=1x_1 =1, xn=g(xn1)=g(g(xn2))=....g(n1)(x1)x_n = g(x_{n-1}) = g(g(x_{n-2}))=.... g^{(n-1)}(x_1) 在本题里面, 对于aa数列,x1=1,g(x)=x+1x2x_1=1, g(x) = x+\frac1{x^2} 对于bb数列,x1=1,g(x)=x(1+1x)3x_1=1, g(x) = x(1+\frac1x)^3 我们可否站在这个层面,即通过分析g(x)g(x)的函数性态,来分析对于给定的点x0x_0g(n)(x0)g^{(n)}(x_0)的稳定性态。 比如,x1=1x_1=1n+n \to +\infty, g(n)(x1)f(n)+cg^{(n)}(x_1) \approx f(n)+c,求f(n)f(n)cc。 更进一步,f(n)f(n)cc是否跟初始值x1x_1无关?

比如,x1=xx_1=xh(x)=g(+)(x)h(x) = g^{(+\infty)}(x) ,求h(x)h(x)

mathe认为:
一般情况是迭代x1=g(x0),..xn=g(xn1)=..=gn(x0)x_1=g(x_0),..x_n=g(x_{n-1})=..=g^n(x_0)有一个不动点x.x^{.},也就是x.=g(x.)x^{.}=g(x^{.}),而且迭代过程会收敛到x.x^{.},我们需要分析收敛情况.
我们可以先做分式线性变换h(x)=1g(1x+x.)x.h(x)=\frac1{g(\frac1x+x^{.})-x^{.}},将不动点变化到无穷远点。而对于本题中,是属于单侧接近x.x^{.},对应可以转化为迭代充分大次数后数值单调增到无穷大 设g(x)=x.+k=1ak(xx.)kg(x)=x^{.}+\sum_{k=1}^{\infty} a_k(x-x^{.})^k, 如果a1<1a_1\lt 1那么会收敛比较快,那么如果其中a1=1a_1=1于是收敛会很慢,其中
h(x)=1k=1akxk=x+k=0bkxkh(x)=\frac 1{\sum_{k=1}^{\infty}\frac{a_k}{x^k}}=x+\sum_{k=0}^{\infty}\frac{b_k}{x^k}
而对于特殊迭代h(x)=x+uxhh(x)=x+\frac u{x^h},可以利用xn+1h+1=xnh+1+uh+x_{n+1}^{h+1}=x_n^{h+1}+uh+\dots来替换成b00b_0\ne 0的模式.
当然如果只分析收敛速度,就不需要展开这么多项,基本只要分析常数项b0b_0的值即可.

Buffalo给出了公式形式:
对于方程an+1=αan+Pk(n)anma_{n+1}=\alpha a_n+\frac{P_k(n)}{a_n^m}Pk(n)=i=0kβiniP_k(n)=\sum_{i=0}^k\beta_i n^i(所有参数都大于零)有以下一般性结论

  1. α>1\alpha\gt 1,则ancαn+...a_n\approx c\alpha^n+...
  2. α=1\alpha=1,则an(m+1k+1βknk+1)1m+1+...a_n\approx (\frac{m+1}{k+1}\beta_k n^{k+1})^{\frac{1}{m+1}}+...
  3. 0<α<10\lt \alpha\lt 1,则an(βknk1α)1m+1+...a_n\approx (\frac{\beta_k n^{k}}{1-\alpha})^{\frac{1}{m+1}}+...

α1\alpha\ge 1时上面的结论没问题。
0<α<10\lt\alpha\lt1时用代换an=(βknk1α)1m+1bna_n=(\frac{\beta_k n^k}{1-\alpha})^{\frac{1}{m+1}}b_n重新定标得到方程bn+1=αbn+1αbnm+o(1n)b_{n+1}=\alpha b_n+\frac{1-\alpha}{b_n^m}+o(\frac{1}{n})
非线性极限普适方程Sn+1=αSn+1αSnmS_{n+1}=\alpha S_n+\frac{1-\alpha}{S_n^m}有不动点,情况比较复杂。随着mm的取值变化数列的“极限”会出现分岔、混沌现象。第一个分岔点在m=1+α1αm=\frac{1+\alpha}{1-\alpha}处。

g(1)(x)=g(x)=(1+x)3x2,g(k+1)(x)=g(g(k)(x))g_{(1)}(x)=g(x)=\frac{(1+x)^3}{x^2}, g_{(k+1)}(x)=g(g_{(k)}(x))我们知道h(g(x))=h(x)+3h(g(x))=h(x)+3 那么我们知道对于所有n>0n>0,方程g(n)(x)=xg_{(n)}(x)=x的解都必然是h(x)h(x)的非解析点. 试做了一下g(5)(x)=xg_{(5)}(x)=x的零点图如下 data
Buffalo做了总结:
做个小结: 考察满足递推关系an+1=αan+Pk(n)anma_{n+1}=\alpha a_n+\frac{P_k(n)}{a_n^m}的初值为xx的数列的极限行为。这里的Pk(n)=i=0kαiniP_k(n)=\sum_{i=0}^{k}\alpha_i n^i,不妨设αk=1\alpha_k=1。所有的参数都是非负正实数,特别是mm也不限于整数。现在根据α\alpha取值范围分情况讨论。

  1. α>1\alpha>1
    此时anh(x)αn1+Qk(n)h(x)mαmn+a_n \approx h(x)\alpha^{n-1}+\frac{Q_k(n)}{h(x)^m}\alpha^{-mn}+\cdots。多项式Qk(n)Q_k(n)满足Qk(n+1)αm+1Qk(n)=α2mPk(n)Q_k(n+1)-\alpha^{m+1}Q_k(n)=\alpha^{2m}P_k(n)。第二项已经很小了,后面的没必要计算下去。
    Pk(n)=1P_k(n)=1,则h(x)h(x)满足h(αx+1xm)=αh(x)h(\alpha x+\frac{1}{x^m})=\alpha h(x),由这个等式可以求得h(x)x(1+αm(αm+11)xm+1i=2m(i1)!R(i,αm+1)j=1i(α(m+1)j1)(αmxm+1)i)h(x)\approx x(1+\frac{\alpha^m}{(\alpha^{m+1}-1)x^{m+1}}-\sum_{i=2}^{\infty}\frac{m}{(i-1)!}\frac{R(i,\alpha^{m+1})}{\prod_{j=1}^{i}(\alpha^{(m+1)j}-1)}(-\frac{\alpha^m}{x^{m+1}})^i)R(i,t)R(i,t)是的tti(i1)21\frac{i(i-1)}{2}-1次多项式,系数除了11次、i(i1)22\frac{i(i-1)}{2}-2次为零之外都是mmi2i-2次正整系数多项式,而且i(i1)21\frac{i(i-1)}{2}-1次系数为(m+i2)!m!\frac{(m+i-2)!}{m!},常数项为j=1i2(im+j)\prod_{j=1}^{i-2}(im+j),二次项系数为(i2)(m+1)j=2i2(im+j)(i-2)(m+1)\prod_{j=2}^{i-2}(im+j)
    前几个R(i,t)R(i,t)R(2,t)=1R(2,t)=1R(3,t)=3m+1+(m+1)t2R(3,t)=3m+1+(m+1)t^2
  2. 0<α<10<\alpha<1
    作代换an=(αknk1α)1m+1bna_n=(\frac{\alpha_k n^k}{1-\alpha})^{\frac{1}{m+1}}b_n可以得到bn+1=αbn+1αbnm+O(1n)b_{n+1}=\alpha b_n+\frac{1-\alpha}{b_n^m}+O(\frac{1}{n})。普适极限方程Sn+1=αSn+1αSnmS_{n+1}=\alpha S_n+\frac{1-\alpha}{S_n^m}有不动点Sn=1S_n=1,随着参数的变化数列的趋势会出现分岔、混沌现象。当0<m<1+α1α0<m<\frac{1+\alpha}{1-\alpha}Sn1S_n \rightarrow 1m>1+α1αm>\frac{1+\alpha}{1-\alpha}时不动点不稳定,开始变成两周期以至更长周期的循环,而且很快就变成混沌的。
  3. α=1\alpha=1
    作代换bn=anm+1b_n=a_n^{m+1}可以得到bn+1=bn+m+1+i=2Cm+1i(Pk(n)bn)i1b_{n+1}=b_n+m+1+\sum_{i=2}^{\infty}C_{m+1}^{i}(\frac{P_k(n)}{b_n})^{i-1},所以bnQk+1(n)+m2lnn+h(Pk,x)+O(lnnn)b_n\approx Q_{k+1}(n)+\frac{m}{2}\ln n+h(P_k,x)+O(\frac{\ln n}{n}),多项式Qk+1(n)Q_{k+1}(n)满足Qk+1(n+1)Qk+1(n)=Pk(n)Q_{k+1}(n+1)-Q_{k+1}(n)=P_{k}(n)。 如果Pk(n)=1P_{k}(n)=1可以进一步得到bn(m+1)n+m2ln(n+m2lnn+h(x)m+1)+h(x)m(2m+1)12((m+1)n+m2lnn+h(x))+m24ln(1+m2lnn+h(x)2(m+1)n)(m+1)x+m2lnn+h(x)7m3+4m248(m+1)2n2+i=3Ri2(m2lnn+h(x))nib_n\approx (m+1)n+\frac{m}{2}\ln(n+\frac{\frac{m}{2}\ln n+h(x)}{m+1})+h(x)-\frac{m(2m+1)}{12((m+1)n+\frac{m}{2}\ln n+h(x))}+\frac{\frac{m^2}{4} \ln(1+\frac{\frac{m}{2}\ln n+h(x)}{2(m+1)n})}{(m+1)x+\frac{m}{2}\ln n+h(x)}-\frac{7m^3+4m^2}{48(m+1)^2 n^2}+\sum_{i=3}^{\infty}\frac{R_{i-2}(\frac{m}{2}\ln n+h(x))}{n^i}