差分方程求解问题

Buffalo

mathe 发表于 2015-12-31 16:14
对于\(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}\fra ...

"Carrier's Rule": Divergent series converge faster than convergent series because they don't have to converge.

Buffalo

wayne 发表于 2015-12-31 16:29
我是用GMP任意精度计算的,程序跑到MinGW64下.
$f[10^8] = 1.135255909724177761096958232161278846371 ...

准确到30位的数值是 1.135255847315503714195394347748
你的数值
1.13525590972417776109695823216127884637104629990206
在第7位就开始出现误差了

wayne

Buffalo 发表于 2015-12-31 16:54
准确到30位的数值是 1.135255847315503714195394347748
你的数值
1.1352559097241777610969582321612 ...

我前面是在计算 \(\Delta(n)=a(n)^3-3n-\ln(n)，\Delta(10^8)\)的精确值，而非收敛值，跟你们的讨论还差一截呢，先算算简单的玩玩。

GMP不适合浮点运算，前面对数是用double初始化的，确实有问题，虽然不会对结果产生累计误差，特改成MPFR，计算得到：
$\Delta(10^8) = 1.1352559097241794977019369379608393295441679846218785948362296416770693539072882376134403370431914736$
$\Delta(10^9) $要计算19分钟：
$\Delta(10^9) = 1.1352558543238998089134184304819402101461919498096029650116492183769582059250023236791929390909381445$
代码如下

1. #include <stdio.h>
   
2. #include <stdlib.h>
   
3. #include <mpfr.h>
   
4. int main(int c,char**v)
   
5. {
   
6.     if(c>1){
   
7.         int n = atoi(v[1]),i;
   
8.         mpfr_set_default_prec(500);
   
9.         mpfr_rnd_t rnd = MPFR_RNDN;
   
10.         mpfr_t an,tmp,tmp2;
    
11.         mpfr_init_set_si(an,1,rnd);
    
12.         mpfr_inits(tmp,tmp2,'\0');
    
13.         for(i=2;i<n+1;++i){
    
14.             mpfr_ui_div(tmp,1,an,rnd);
    
15.             mpfr_pow_ui(tmp,tmp,2,rnd);
    
16.             mpfr_add(an,an,tmp,rnd);
    
17. //            mpfr_printf("f[%d] = %.100RNf\n",i,an);
    
18.         }
    
19. 
    
20.         mpfr_pow_ui(tmp,an,3,rnd);
    
21.         mpfr_sub_ui(tmp,tmp,3*n,rnd);
    
22.         mpfr_set_ui(tmp2,n,rnd);
    
23.         mpfr_log(tmp2,tmp2,rnd);
    
24.         mpfr_sub(tmp,tmp,tmp2,rnd);
    
25.         mpfr_printf("f[%d] = %.100RNf\n",n,tmp);
    
26.         printf("This is MPFR %s!\n",mpfr_version);
    
27.         mpfr_clears(an,tmp,tmp2,'\0');
    
28.     }
    
29.     return 0;
    
30. }
    

复制代码

Buffalo

绕开数列通式的渐近式直接快速求极限值的渐近展开式方法。
根据前面的分析，我们知道$h(b_i)=h(b_1)+3i-3$，因为$b_2=b_1+3+\frac{3}{b_1}+\frac{1}{b_1^2}$，所以$h(x)$应该满足函数方程$h(x+3+\frac{3}{x}+\frac{1}{x^2})=h(x)+3$，从这个函数方程即可快速求得极限值$a(x)$的渐近式：前面的几项为
$h(x)~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)!*3^n$化为整数实在是看不出规律，连符号都无法预测。

mathe

符号比较错乱，尽力统一一下吧。
\(a(n+1)=a(n)+\frac{1}{a(n)^2},b(n)=a(n)^3,\Delta(n)=b(n)-3n-\ln(n)\)
另外，为了便于区分在\(a(1)=x\)时，我们可以将\(\Delta(n)\)写成\(\Delta_x(n)\),类似有\(a_x(n),b_x(n)\)
定义\(h(x)=h_x=\lim_{n->\infty}\Delta_x(n)\)

mathe

$u_m(n)=b_n-3n-\ln(n)-a-\sum_{k=1}^m\frac{P_k(a+\ln(n))}{n^k}$
如果能够证明$\lim_{n->\infty}u_m(n)n^m=0$那么就可以证明这个式子是$b_n$的渐进式，其中$a=\lim_{n->\infty}b_n-3n-\ln(n)$

wayne

换种表达应该更简洁吧，设$c_n = 1/a_n^3 = 1/b_n$，那么 $c_{n} = c_{n-1}/(1+c_{n-1})^3$

mathe

记\(f(x)=x+3+\frac{3}{x}+\frac{1}{x^2},f_1(x)=f(x),f_{k+1}(x)=f(f_k(x))\)。那么\(f_k(x)\)相当于保持数列\(b(n)\)递推关系不变，但是第零项\(b(0)\)改变成$x$,那么第$k$项就变为\(f_k(x)\)了。
设\(y=g(a,m)\)为使得\(\displaystyle\lim_{n->\infty}f_n(y)-3(n+m)-\ln(n+m)=a\)的数而且\(y\ge 2\) 。
根据前面的分析过程，不管数列\(f_n(y)\)中第零项\(y\)是多少，极限\(\displaystyle\lim_{n->\infty}f_n(y)-3n-\ln(n)\)都是存在的，
于是对于给定的m极限\(\displaystyle\lim_{n->\infty}f_n(y)-3(n+m)-\ln(n+m)\)也存在只是和前一个极限差了常数\(3m\)。
于是\(y=g(a,m)\)表示如果保持数列\(b(n)\)的递推关系，其中第\(m\)项\(b(m)\)为\(y\),那么\(\displaystyle\lim_{n->\infty}b(n)-3n-\ln(n)\)的极限正好为\(a\)。
而比较有意思的是这时候的\(m\)我们可以推广到任意一个实数而不仅仅是整数，从而将数列\(b_a(m)\)推广为一个定义在正实数范围以\(m\)为自变量的函数
于是容易得出\(g(a,m)=g(a+3d,m-d)\),而且对于我们以前定义的\(b_a(m)=g(a,m)\)
而根据前面分析应该有$g(a,n)~3n+\ln(n)+a+\sum_{k=1}^{\infty}\frac{P_k(\ln(n)+a)}{n^k}$
由于\(g(a,n)=g(a-3d,n+d)\),如果我们选择适当d使得\(\ln(n+d)+a-3d=0\),记其中\(m=n+d\),得到
那么结果可以简化为$g(a,n)=g(a-3d,m) "~" 3m+\sum_{k=1}^{\infty}\frac{P_k(0)}{m^k}$
由此，对于任意$n,b_a(n)=g(a,n)$, 如果我们可以先利用上面公式反解出对应的m，就可以求出对应的$d=m-n$,然后利用$a=3d-log(m)$得出$a$
于是问题就变成求渐进式$s(x)~3x+\sum_{k=1}^{\infty}\frac{P_k(0)}{x^k}$,以及其逆函数
然后我们假设当将n替换成n+1时，对应的$m$变化为$m'$，于是根据m的定义必然有$m'-\frac{\ln(m')}{3}=m+1-\frac{\ln(m)}{3}$
这里$m'$可以看成是$m$的隐函数，显然可以看出在$m -> \infty$时,$m' = m+1+o(1)$,代入递推式$m'=m+1+\frac{\ln(\frac{m'}{m})}{3}$,可以看出对于一般情况的f,这里$m',m$之间关系式直同f中常数项相关
得出$m' = m+1+\frac{1}{3m}+o(\frac{1}{m})$,反复迭代可以得出$m'=J(m)=m+1+\frac{1}{3m}-\frac{1}{18m^2}-\frac{1}{54m^3}+\frac{7}{324m^4}+...$
系数依次为$1/3,-1/18,-1/54,7/324,-31/4860,-107/29160,2833/612360,-621/459270,...$
然后再次利用公式$b_a(n+1)=b_a(n)+3+3/{b_a(n)}+1/{b_a(n)^2}$,并且$b_a(n+1)=s(m'),b_a(n)=s(m)$,将$m'=J(m)$代入，
也就是求$s(J(m))=s(m)+3+3/{s(m)}+1/{s(m)^2}$
比较系数，就可以得出${P_k(0)}$的值,
分别为$-5/18,-1/6,-157/1458, -13327/174960,-444191/7873200,-14533499/330674400,-776933609/20832487200$
从前面的结果来看都是负数而且绝对值小于1，也就是级数很可能是收敛的
于是$H(x)=\frac{1}{s(1/x)}$满足H(0)=0而且展开式非常容易算出，于是H的逆函数在0的展开式也可以计算出来，由此可以得出s的逆函数的展开式，也是罗兰级数形式
https://en.m.wikipedia.org/wiki/Lagrange_inversion_theorem
可以用来计算更高精度结构，然后结果再倒数一下可以得到s的逆函数形式，也就是可以用$b_a(n)$直接计算对应的$m$,最后可以推导出对应的$a$
最后m的公式应该是$1/3*x+5/(18*x)+1/(2*x^2)+239/(324*x^3)+5227/(6480*x^4)+...$

wayne

设$c_n = 1/a_n^3 = 1/b_n$，那么 $c_{n} = c_{n-1}/(1+c_{n-1})^3$ ，为了方便分析的引入，可以进行适当延拓。
设$g(x+t) ={g(x)g(t)}/(1+g(x t))^3$,其中$g(1)=1$，那么 $c(n)$是 $g(x)$的一个特解。
于是标记$g(0)=a, g'(0)=b$, 对$t$进行求导，得到微分方程：$g'(x) = ((1+ a - 3 a x) b g(x))/(1 + a)^4$
解得 $g(x) =a exp(\frac{b (-{3 a x^2}/2+a x+x)}{(a+1)^4})$ ，  其中 $a$是 $(1+a)^3=a$的根，根据$g(1)=1$得知， $b$是 $b(1-a/2) + (1+a)^4ln(a)=0$的根

=================
有点不对劲，。。。。

Buffalo

上面给出了初值$x$很大时数列极限$H(x)$的渐近式$H_m(x)$，如果$x$非常小，则下一项$\frac{(1+x)^3}{x^2}$就会很大，用这个很大的值就可以算出很小的初值数列极限$h(x)$的渐近式$h_{2m}(x)=H_m(\frac{(1+x)^3}{x^2})-3$。$h_{2m}(x)$比$H_m(x)$好看多了，至少系数符号看起来是严格交错的，而且乘以$3*n!$之后是整数。
现在的问题是：如何有效整合这两个渐近式，给出一个有效计算初值在$1$附近的数列极限的公式？