数学研发论坛

 找回密码
 欢迎注册
查看: 39563|回复: 81

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

[复制链接]
发表于 2015-12-29 18:33:35 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
对于\(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}\)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-29 19:08:24 | 显示全部楼层
\(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}\)

点评

得到的结果还需要进一步精密估计,数值分析你的结果不如1#彩色の夢∩o∩)精确  发表于 2015-12-29 21:12

评分

参与人数 1威望 +6 金币 +4 收起 理由
数学星空 + 6 + 4 赞一个!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-29 20:10:19 | 显示全部楼层
然后我们可以利用积分估算得出
\(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)\)只略微大了一些
同样要数值计算这个极限值,也可以先计算前面若干项,然后对于余下部分求和式和近似公式求差值,而差值部分误差可以非常容易估计,这样就可以得出极限值一个较好的估计范围了。

评分

参与人数 1金币 +12 贡献 +12 经验 +12 收起 理由
数学星空 + 12 + 12 + 12 很给力!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-30 14:46:54 | 显示全部楼层
本帖最后由 Buffalo 于 2015-12-30 17:05 编辑

首先令$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},找不到规律,暂停。

点评

不太理解每一项为什么都含3n+ln(n)+a?  发表于 2015-12-30 21:30
有些矛盾在里面  发表于 2015-12-30 21:27

评分

参与人数 1金币 +12 贡献 +12 经验 +12 鲜花 +12 收起 理由
数学星空 + 12 + 12 + 12 + 12 神马都是浮云

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-31 12:14:09 | 显示全部楼层
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})−\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}$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-31 14:24:08 | 显示全部楼层
充分利用。
假设已经求得初值为$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}$,应该去掉。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-31 16:14:29 | 显示全部楼层
对于\(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\)

点评

嗯!多谢指点  发表于 2016-1-1 08:14
y->y/(1+y)就是n->n+1 X->X+ln(1+y)就是a+ln(n)->a+ln(n+1)=a+ln(n)+ln(1+1/n)  发表于 2015-12-31 21:57
代换y->y/(1+y),X->X+ln(1+y),很神奇,是如何得到此变换的啊?  发表于 2015-12-31 21:43

评分

参与人数 1威望 +12 金币 +20 贡献 +12 收起 理由
数学星空 + 12 + 20 + 12 神马都是浮云

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-31 16:17:17 | 显示全部楼层
分析能力没大家强,我就怒算至$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!

点评

额,我是算b(10^8) 的精确值,而非你们的收敛值。。。  发表于 2015-12-31 23:25
嗯嗯,只是 tail了一下最后的10个,肉眼观察,随口一说,没细细追究  发表于 2015-12-31 18:10
只是变化很慢,不代表有这么高精度  发表于 2015-12-31 18:06
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-31 16:24:47 | 显示全部楼层
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}$,直接算效率是很低的,而且算这么多次舍入误差累积起来不得了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-12-31 16:29:01 | 显示全部楼层
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. //        printf("This is GMP %s!\n",gmp_version);
  27.         mpf_clears(bn1,bn2,tmp,tmp2,'\0');
  28.     }
  29.     return 0;
  30. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2022-7-1 10:02 , Processed in 0.221166 second(s), 18 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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