找回密码
 欢迎注册
楼主: 数学星空

[讨论] 差分方程求解问题

[复制链接]
发表于 2015-12-31 16:41:45 | 显示全部楼层
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.

点评

主要问题是在于如何说明这中方法得出的结果是可靠的。  发表于 2015-12-31 17:27
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-31 16:54:24 | 显示全部楼层
wayne 发表于 2015-12-31 16:29
我是用GMP任意精度计算的,程序跑到MinGW64下.
$f[10^8] = 1.135255909724177761096958232161278846371 ...


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

点评

10^9跑了19分钟  发表于 2016-1-1 09:21
我是昨天刚接触这个题目,正好想练练MSYS2/MinGW64 + Qt 64位下编译64位程序的效果。所以没看你们的讨论,直接计算数组的精确值。。。  发表于 2016-1-1 09:21
1.13525584731550371419539434774792944 原来是wayne的数据算错了,利用他的结果可以得出这个值  发表于 2016-1-1 07:46
我用四次逼近的公式逆求a好像精度也不高  发表于 2015-12-31 21:56
1.1352558473155019775904214306  发表于 2015-12-31 21:39
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-31 23:02:49 | 显示全部楼层
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.         mpfr_pow_ui(tmp,an,3,rnd);
  20.         mpfr_sub_ui(tmp,tmp,3*n,rnd);
  21.         mpfr_set_ui(tmp2,n,rnd);
  22.         mpfr_log(tmp2,tmp2,rnd);
  23.         mpfr_sub(tmp,tmp,tmp2,rnd);
  24.         mpfr_printf("f[%d] = %.100RNf\n",n,tmp);
  25.         printf("This is MPFR %s!\n",mpfr_version);
  26.         mpfr_clears(an,tmp,tmp2,'\0');
  27.     }
  28.     return 0;
  29. }
复制代码

点评

Buffalo已经在9楼做了误差分析,我不是在计算那个收敛值,只是计算精确值,没别的意思  发表于 2016-1-1 09:33
死算\delta(10^n),精确到小数点后10位看来都很困难.  发表于 2016-1-1 09:25
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-1-1 06:35:21 | 显示全部楼层
绕开数列通式的渐近式直接快速求极限值的渐近展开式方法。
根据前面的分析,我们知道$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$化为整数实在是看不出规律,连符号都无法预测。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-1-1 08:01:47 | 显示全部楼层
符号比较错乱,尽力统一一下吧。
\(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)\)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-1-1 09:15:35 | 显示全部楼层
$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)$

点评

Buffalo 应该是靠展开式和误差来计算的  发表于 2016-1-1 09:36
mathe利用是$\Delta$的精确值反推的  发表于 2016-1-1 09:35
如何验算a的精确值,Buffalo 给出 1.1352558473155037141...与你给出的1.1352558473155019775904214306不一致啊?  发表于 2016-1-1 09:27
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-1-1 09:50:20 | 显示全部楼层
换种表达应该更简洁吧,设$c_n = 1/a_n^3 = 1/b_n$,那么 $c_{n} = c_{n-1}/(1+c_{n-1})^3$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-1-1 21:32:01 | 显示全部楼层
记\(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)+...$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-1-1 23:05:06 | 显示全部楼层
设$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$的根

=================
有点不对劲,。。。。

点评

选择双变量方程是有规可循的,比如对称性,初始值。。。, 没选好是回不去的。  发表于 2016-1-3 11:03
比如: f(x+y) =f(x)*f(y), f(xy) = f(x)+f(y) ....  发表于 2016-1-3 11:01
确实是必要非充分,但很多的双变量的单元函数方程都可以通过微分手段得到答案  发表于 2016-1-3 11:00
不对劲的原因在于后面的单变量微分方程只是双变量方程$g(x+t)=\frac{g(x)g(t)}{(1+g(xt))^3}$的必要非充分条件,弱化了很多,回不去也是理所当然的。不管怎么选择双变量方程都基本不可能走通这条路。  发表于 2016-1-2 04:25
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-1-2 10:56:37 | 显示全部楼层
本帖最后由 Buffalo 于 2016-1-2 11:20 编辑

上面给出了初值$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$附近的数列极限的公式?

点评

最小初值是2,由于函数$x+3+3/x+1/x^2$在$x<2$式减,所以只要迭代一次反尔会有更大的值,所以只要考虑$x>=2$区间即可  发表于 2016-1-2 21:59
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-28 21:10 , Processed in 0.054272 second(s), 16 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表