差分方程求解问题

数学星空

对于\(a(n+1)=a(n)+\frac{1}{a(n)^2},a(1)=1\),

我们可以利用微分方程求解\(a(n)\)的第一阶项.

设\(a(n)=y\)

又\(a(n+1)-a(n)=a'(n)\)

则有\(y'y^2=1\)

求解得到:\(a(n)=y \approx (3n)^{\frac{1}{3}}\)

接下来,我们可设

\((a(n))^3=3n+a+b\ln(n)+\frac{c}{n}+\frac{d}{n^2}+\frac{e}{n^3}+\frac{f}{n^4}+...\)

本主题的任务是如何求解\(a,b,c,d,e\)?

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

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

\((3n)^{1/3}<a(n)<(3n)^{1/3}+(3n)^{-1/3}\)

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

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

mathe

\(a(n+1)=a(n)+\frac{1}{a(n)^2}\)
得出
\(a(n+1)^3=a(n)^3+3+\frac{3}{a(n)^3}+\frac{1}{a(n)^6}\)
归纳容易得出对于\(n\ge 2\)有\(a(n)^3\ge 3n+2\)于是重新代入上面得出
\(a(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}\)

mathe

然后我们可以利用积分估算得出
\(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)^3 \ge 8+3(n-2) +3\sum_{k=2}^{n-1} \frac{1}{3n+2+\ln(n)}\)
而且容易看出后面求和式\(3\sum_{k=2}^{n-1}\frac{1}{3n+2}\)的差在\(n\)趋向无穷收敛到一个常数
所以可以知道存在常数\(C\)使得\(a(n)^3 \ge 3n+\ln(n)+C\)
也就是我们可以非常轻松证明\(a(n)^3-3n-\ln(n)\)有界，进一步分析还可以得出其极限存在，但是求这个极限比较有难度
而数值计算\(\Delta(n)=a(n)^3-3n-\ln(n)\),有
$\Delta(10)  =1.220832$
$\Delta(100)=1.151530$
$\Delta(1000)=1.137657$
$\Delta(10000)= 1.135573$
$\Delta(100000)=1.135295$
$\Delta(1000000)= 1.135261$
$\Delta(10000000)=1.135256$
所以可以看出1＃彩色の夢∩o∩)估计公式中常数项偏离比较大，倒是本楼开头估计式中对应常数项\(\ln(3)+\frac{1}{8}+2-\ln(8)\)只略微大了一些
同样要数值计算这个极限值，也可以先计算前面若干项，然后对于余下部分求和式和近似公式求差值，而差值部分误差非常容易估计，这样就可以得出一个较好的估计范围了。

Buffalo

首先令$b_n=a_n^3$，得到$b_{n+1}=b_n+3+\frac{3}{b_n}+\frac{1}{b_n^2}$。渐近式前三项当然是$3n+\ln n+a$，但是后面的项可不是$\frac{1}{n}$的简单级数，而应该是$\sum_{i=1}^{\infty}\frac{P_i(ln n,a)}{n^i}$，这里的$P_i(\ln n,a)$是$\ln n$的$i$次多项式（事实上同时也是$a$的$i$次多项式）。用待定系数法，令$f_m(n)=3n+\ln n+a+\sum_{i=1}^{\m}\frac{P_i(ln n,a)}{n^i}$让mathematica在$n=\infty$处展开$f_{m}(n+1)-f_m(n)-3-\frac{3}{f_m(n)}-\frac{1}{f_m(n)^2}$到$\O(\frac{1}{n^{m+1}})$逐次得到$P_i(\ln n,a)$。计算发现$P_i(\ln n,a)$实际上很神奇地仅仅是是$\ln n+a$的$i$次多项式，而且系数有很明显的规律，可以把一部分项合并降低次数，目前做到$b_n~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}，找不到规律，暂停。

Buffalo

Buffalo 发表于 2015-12-30 14:46
首先令$b_n=a_n^3$，得到$b_{n+1}=b_n+3+\frac{3}{b_n}+\frac{1}{b_n^2}$。渐近式前三项当然是$3n+\ln n+a$ ...

还是找不到规律。不过目前的结果用于快速计算任意初值的数列极限值已经足够了：已经知道$b_n~3n+\lnn+a+ln(1+\frac{\ln n+a}{3n})&#8722;\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}}$，从初值出发计算$N$项，定义$\Delta(n)=b_n-3n-\ln n$，$\Delta(N)$就是极限a的粗略值$d_0$，然后令$d_{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(\frac{\ln N}{N^3})$。
费点功夫把展开式多计算几项就可以得到更好的逼近效果$\O(\frac{(\ln N)^i}{N^{i+2}})$。
作为例子，选初值为1，计算100项，迭代几次后得到极限值为1.1352567349497837，计算10000项，迭代几次后得到极限值为1.1352558474127066，两者相差$8.9*10^{-7}$，小于估计的$\O(\frac{\ln N}{N^3})~5*10^{-6}$

Buffalo

充分利用。
假设已经求得初值为$b_1$的数列的通式$b_n=f(n,b_1)~3n+\ln n+h_{b_1}+R(n,a(b_1))$，现在来看对极限值$h_{b_1}$我们能了解到什么程度。
递推方程是平移不变的，因此$f(n,b_i)=f(n,f(i,b_1))=f(n+i-1,b_1)$，所以有$h_{b_i}=\lim_{n->\infty}(f(n,b_i)-3n-\ln n)=\lim_{n->\infty}f(n+i-1,b_1)-3n-\ln n=h_{b_1}+3i-3$，这是个严格的等式。由$b_i=f(i,b_1)~3i+\ln i+h_{b_1}+R(i,h_{b_1})$可以迭代求得$3i~b_i-\ln\frac{b_i}{3}-h_{b_1}+...$，最后得到$h_{x}~x-\ln \frac{x}{3}-3-\frac{\ln (x/3)}{x}+\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)$里含有巨量的$\ln n+a$，可是后面的项里面不仅$h_{b_1}$合情合理地消失了，连$\ln x$也无影无踪。
对于很大的初值$b_1$，可以直接用这个表达式作为近似值。如果$b_1$不太大，可以先迭代$i-1$次得到较大的$b_i$，算出$h_{b_i}$，然后用等式$h_{b_1}=h_{b_i}-3i+3$计算。

补充内容 (2016-1-1 06:11):
$h_{x}$中不知道怎么多写了一项$-\frac{\ln \frac{x}{3}}{x}$，应该去掉。

mathe

对于\(b_{n+1}=b_n+3+\frac{3}{b_n}+\frac{1}{b_n^2}\),假设
\(b_n=3n+\ln(n)+a+\sum_{k=1}^{+\infty}\frac{P_k(\ln(n)+a)}{n^k}\)

我们先查看\(\frac{P_k(\ln(n+1)+a)}{(n+1)^k}=\frac{P_k(\ln(n)+a+\ln(1+\frac{1}{n}))}{n^k*(1+\frac{1}{n})^k}\)
这个式子展开后是一个无穷级数，可以写成形如\(\sum_{h=k}^{+\infty}\frac{Q_{k,h}(\ln(n)+a)}{n^h}\)的级数，其中
每个\(Q_{k,h}(.)\)都是次数不超过k的多项式,而且同\(a\)无关。
所以从形式上，这种形式的级数应该是可以用于这种递推式的，但是得出的结果不一定是收敛的，如同\(\Gamma\)函数的Stirling级数就是不收敛的，而只是一种渐近式。
为了计算方便，我们可以将\(b_n=3n+\ln(n)+a+\sum_{k=1}^{+\infty}\frac{P_k(\ln(n)+a)}{n^k}\)中n用$\frac{1}{y}$代替，而\(\ln(n)+a\)用\(X\)代替，写成形式级数
\(b_n=\frac{3}{y}+X+\sum_{k=1}^{+\infty}P_k(X)y^k\)
而其中如果我们将\(y\)替换为\(\frac{y}{1+y}\),然后再将\(X\)替换为\(X+\ln(1+y)\),就可以得出\(b_{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 > 3/y+X+(b00+b01*X)*y+(b10+b11*X+b12*X^2)*y^2+O(y^3)
%1 = 3*y^-1 + X + (b01*X + b00)*y + (b12*X^2 + b11*X + b10)*y^2 + O(y^3) //这个即$b_n$二阶近似
(16:49) gp > subst(%1,y, y/(1+y))
%3 = 3*y^-1 + (X + 3) + (b01*X + b00)*y + (b12*X^2 + (-b01 + b11)*X + (-b00 + b10))*y^2 + O(y^3)
(16:51) gp > subst(%3,X,X+log(1+y))
%4 = 3*y^-1 + (X + 3) + (b01*X + (b00 + 1))*y + (b12*X^2 + (-b01 + b11)*X + (-b00 + (b01 + (b10 - 1/2))))*y^2 + O(y^3)//这个为$b_{n+1}$的二阶近似
(16:52) gp > %4-%1-3-3/%1-1/%1^2
%6 = ((-b01 + 1/3)*X + (-b00 + (b01 - 11/18)))*y^2 + O(y^3)
由此得出b01=1/3,b00=-5/18
类似可以得出
\(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}+...\)

公式:
\(\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}\)
\((\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\)

wayne

分析能力没大家强,我就怒算至$10^7$,发现稳定的数据是  $1.1352563946493.......$,mathe的数据是$\Delta(10000000)=1.135256$,
f[9999992] = 1.13525639464972320465721775596975762052600000000000
f[9999993] = 1.13525639464967102900607028356514436617500000000000
f[9999994] = 1.13525639464961951149190589004304185411700000000000
f[9999995] = 1.13525639464956865211672457638613588410000000000000
f[9999996] = 1.13525639464951489816884754247600004090200000000000
f[9999997] = 1.13525639464946535507763238909718281733600000000000
f[9999998] = 1.13525639464941291741772151443072277723700000000000
f[9999999] = 1.13525639464936113790479371755988483819800000000000
f[10000000] = 1.13525639464931001654084899646744048433600000000000
Hello World,6.0.0!

Buffalo

wayne 发表于 2015-12-31 16:17
分析能力没大家强,我就怒算至$10^7$,发现稳定的数据是  $1.1352563946493.......$,mathe的数据是b(10000000 ...

注意到$b_n~3n+\ln n+a+\O(\frac{\ln n}{n})$，算到$10^7$项也只能把误差控制在$10^{-6}$，直接算效率是很低的，而且算这么多次舍入误差累积起来不得了。

wayne

Buffalo 发表于 2015-12-31 16:24
注意到$b_n~3n+\ln n+a+\O(\frac{\ln n}{n})$，算到$10^7$项也只能把误差控制在$10^{-6}$，直接算效率是 ...

我是用GMP任意精度计算的,程序跑到MinGW64下.
$f[10^8] = 1.13525590972417776109695823216127884637104629990206$

1. #include <stdio.h>
   
2. #include <gmp.h>
   
3. #include <math.h>
   
4. #include <stdlib.h>
   
5. int main(int c,char**v)
   
6. {
   
7.     if(c>1){
   
8.         int n = atoi(v[1]);
   
9.         mpf_set_default_prec(200);
   
10.         mpf_t bn1,bn2,tmp,tmp2;
    
11.         mpf_init_set_si(bn1,1);
    
12.         mpf_inits(bn2,tmp,tmp2,'\0');
    
13.         for(int i=2;i<n+1;++i){
    
14.             mpf_ui_div(tmp,1,bn1);
    
15.             mpf_add_ui(tmp,tmp,1);
    
16.             mpf_pow_ui(tmp,tmp,3);
    
17.             mpf_mul(tmp,bn1,tmp);
    
18.             mpf_set(bn2,bn1);
    
19.             mpf_set(bn1,tmp);
    
20.             mpf_sub_ui(tmp,tmp,3*i);
    
21.             mpf_set_d(tmp2,log(i));
    
22.             mpf_sub(tmp,tmp,tmp2);
    
23. //            gmp_printf("f[%d%] = %.50Ff\n",i,tmp);
    
24.         }
    
25.         gmp_printf("f[%d%] = %.50Ff\n",n,tmp);
    
26. 
    
27. //        printf("This is GMP %s!\n",gmp_version);
    
28.         mpf_clears(bn1,bn2,tmp,tmp2,'\0');
    
29.     }
    
30.     return 0;
    
31. }
    

复制代码