TSC999 发表于 2022-6-15 18:18:30

数列 a(1)=1, a(n+1)=a(n)[1-a(n)/3], 求证 5/2<100a(100)<3

问题 ①: 2022 年浙江高考数学卷第 10 题:



问题 ②: 如何用软件 mathematica 在短时间内求出该数列第 100 项的近似数值? 有一个指令可以使用,但是时间太长,例如计算第 25 项,需要用 18 秒。
想算到第 100 项可能需要几十个小时甚至几个月?

Timing[
N@RecurrenceTable[{a == a (1 - a/3), a == 1},
   a, {n, 25, 25}]
]

计算结果是 a(25)≈0.101709。

还有一个程序计算稍快一点,计算第 25 项用时 12 秒,虽有改进,也不咋的:

Timing = 1; a = 2/3;
a := a = a - a a/3;
N@a
]

上面这两个程序能不能改进到用时大量缩短,达到计算第 100 项不超过 10 分钟?

wayne 发表于 2022-6-15 19:59:52

数学归纳法证明$a_n<=\frac{3}{n+2}$.对于 $n=1$成立。 对于一般情况,因为$f(x) =x -1/3x^2$在$(0,3/2)$内单调递增。 所以如果 $a_n< \frac{3}{n+2}$,那么,$a_{n+1}=f(a_n)<f(\frac{3}{n+2}) <3/{n+3}$ 确实恒成立,所以,对于任意的$n>1$,$a_n<\frac{3}{n+2}$

同样,可以证明$na_n$单调递增。因为$(n+1)a_{n+1}-na_n =a_n(1-\frac{n+1}{3}a_n)>=\frac{1}{n+2}a_n >0$

With[{m = 5000}, RecurrenceTable[{a == a (1 - a/3), a == 1`100}, a, {n, m - 10, m}]] // Reverse

更进一步的,得到$a_n <=\frac{3}{n+2}-\frac{275 (n-1)}{54 (n+2)^3}$ 等号在$n=1,3$的时候取得。$n=12$取得最大残差。

lsr314 发表于 2022-6-15 22:29:12

需要证明$na_n$是递增的

mathe 发表于 2022-6-16 14:41:38

查看:
https://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=8790&pid=61918
先定义$b_n=\frac1{a_n}$,
记\(h(x)=\frac{x^2}{x-\frac13}\)得到\(b_{n+1}=h(b_n)\)

wayne 发表于 2022-6-16 17:33:13

受楼上启发,设$c_n=b_n-1/3 = 1/a_n-1/3$,那么,$c_1=2/3, c_n =\frac{n+1}{3} +1/9\sum_{k=1}^{n-1}1/c_k$

TSC999 发表于 2022-6-16 23:59:15

mathe 发表于 2022-6-18 08:38:12

https://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=8790&pid=92433&fromuid=20
$v(x)=x^2/(x-1/3)=x+1/3+1/{9x}+1/{27x^2}+1/{81x^3}+\cdots$
于是我们通过$a_1=1, a_{n+1}=v(a_n)$递推式得出来的数列和本题中数列互为倒数。
对比链接中记号有$b_0=1/3,b_1=1/9$

于是$1/3 \ln(m)-1/3 m = 1/3 \ln(m')-1/3 (m'-1)$
即$\ln(m')-m'+1=\ln(m)-m$, 即$ m' =m +1 +\ln({m'}/m)$
于是$m'=J(m)=m+1 + \ln(\frac{J(m)}m)$
使用迭代计算
? x+O(x^2)
%8 = x + O(x^2)
? log(1+x+x*%)
%9 = x + 1/2*x^2 + O(x^3)
? log(1+x+x*%)
%10 = x + 1/2*x^2 - 1/6*x^3 + O(x^4)
? log(1+x+x*%)
%11 = x + 1/2*x^2 - 1/6*x^3 - 5/12*x^4 + O(x^5)
? log(1+x+x*%)
%12 = x + 1/2*x^2 - 1/6*x^3 - 5/12*x^4 - 1/20*x^5 + O(x^6)
? log(1+x+x*%)
%13 = x + 1/2*x^2 - 1/6*x^3 - 5/12*x^4 - 1/20*x^5 + 49/120*x^6 + O(x^7)
? log(1+x+x*%)
%14 = x + 1/2*x^2 - 1/6*x^3 - 5/12*x^4 - 1/20*x^5 + 49/120*x^6 + 15/56*x^7 + O(x^8)
? log(1+x+x*%)
%15 = x + 1/2*x^2 - 1/6*x^3 - 5/12*x^4 - 1/20*x^5 + 49/120*x^6 + 15/56*x^7 - 457/1260*x^8 + O(x^9)
? log(1+x+x*%)
%16 = x + 1/2*x^2 - 1/6*x^3 - 5/12*x^4 - 1/20*x^5 + 49/120*x^6 + 15/56*x^7 - 457/1260*x^8 - 49/90*x^9 + O(x^10)
? log(1+x+x*%)
%17 = x + 1/2*x^2 - 1/6*x^3 - 5/12*x^4 - 1/20*x^5 + 49/120*x^6 + 15/56*x^7 - 457/1260*x^8 - 49/90*x^9 + 391/2016*x^10 + O(x^11)
得到$J(1/x)=1/x+1+x + 1/2*x^2 - 1/6*x^3 - 5/12*x^4 - 1/20*x^5 + 49/120*x^6 + 15/56*x^7 - 457/1260*x^8 - 49/90*x^9 + 391/2016*x^10 + O(x^11)$
然后得到渐进式$s(x)\approx 1/{3x}-(x/6+{5x^2}/18+{13x^3}/27+{2009x^4}/2160+{6973x^5}/3600+{1923143x^6}/453600+{30503779x^7}/3175200+{850720483x^8}/38102400+{1209787591x^9}/22861440+{3204878097563x^10}/25147584000+{43027959566971x^11}/138311712000+\cdots)$

\(s(\frac1{J(m)}) = v(s(\frac1m))\), 而且\(s(x)=1/{3x}+P_1 x+P_2x^2+...\)
根据现有展开结果,我们可以先计算$a_{1000}=0.0029760594191188835248176874067736714838$, 计算\(\frac1m=s^{-1}(1/{a_{1000}})=0.00099201931744002304343664937365698435425\)。
所以对应\(m=1008.0448862433159869955\), 那么对应极限应该约等于$c_{a_0} =0.37637275522427338461471882180513818568$
但是代回n=2000情况,需要计算对应的m为$1/3 \ln(m)+c_{a_0}-1/3(m-n)=0$,对应m=2008.7343784059206482994905154424269013
然后计算对应数列倒数值为\(s(1/m)=669.57804309542223061154167143880331827\),得出原始数列第2000项为\(1/{s(1/m)}=0.0014934778855307969063639720682167818331\).
直接迭代计算的结果为0.0014934778855307969063639720682167818339,差别应该是舍入误差引起的。

于是对于本题中数列中较大的n,可以通过如下方法计算第n项:
i) 解关于m的方程 $\ln(m)-(m-n)+3*0.37637275522427338461471882180513818568=0$,得到实数m (m略大于n,可以在(n,2n)中使用二分法找m)
ii) 计算$1/{s(1/m)}$, 即为原数列第n项。

mathe 发表于 2022-6-19 06:56:55

如果我们尝试把上面的s(x)的渐进展开式试着转化为连分式形式,会得到
$s(x)=1/{3x}+[-x/6, -{5x}/3, -x/15, -{619x}/210, {419x}/120,{577191x}/9077635, - {4163475702921x}/1746509741095,
{39018927590570016701x}/84109224655563621885, {8320244433144139792056998559830x}/8142103810321752087542998506339,
- {146451852647634675917906267906355624753965x}/57252060088539150878523070380820686585512]$
比较有意思的是不像前面渐进展开式,系数绝对值越来越大,连分式展开式中每项的值都不是特别大。

mathe 发表于 2022-6-19 10:58:46

47#另外一种推导方法有\(C(x)\approx 3x-\ln(x)+\alpha_0 + \sum_{h=1}^{\infty} \frac{\alpha_h}{x^h}\), 而且$C(v(x))=C(x)+1$
设$x=1/t, C(x)=3x-\ln(x)+\alpha_0 +h(\frac1x)$,得到\(h(t)-h(t-\frac{t^2}3)=\frac t{3-t}-\ln(\frac3{3-t})\)
由此得到\(h(t)=\frac t6+\frac{t^2}{27}+\frac{13t^3}{972}+\frac{113t^4}{19440}+\frac{1187t^5}{437400}+\frac{877t^6}{688905}+\frac{14569t^7}{25719120}+\frac{176017t^8}{793618560}+\frac{1745717t^9}{26784626400}+\frac{88217t^{10}}{15345358875}-\frac{147635381t^{11}}{19445638766400}-\frac{3238110769t^{12}}{827323540243200}+\frac{63045343657t^{13}}{37643221081065600}+\frac{24855467017t^{14}}{7125970336860375}+O(t^{15})\),看上去这是一个收敛的级数。
而我们取倒数以后数列第100项为:35.261925907006479247389066690274451303,
计算$C(35.261925907006479247389066690274451303)=102.22773055434092984523940170232737503+\alpha_0$,计算中迭代式最后一项值为8E-28,说明至少精确到小数点后27位。所以我们得到$C(a_{1000})=900+102.22773055434092984523940170232737503+\alpha_0$,二分法(牛顿法也应该可行)求解得到$a_{1000}=336.01479647071971871965605047571202386$,
如果回到原始数列还需要求一次倒数,就是0.0029760594191188835248176874067737027128精度超过了27位。
可以看出上面计算过程中$\alpha_0$的真正取值我们不需要关心。

现在我们再查看\(h(t)-h(t-\frac{t^2}3)=\frac t{3-t}-\ln(\frac3{3-t})\),设\(h(t)=\sum_{k=1}^{\infty} h_k t^k\)
于是我们有
\(h_k - \sum_{s=0}^{\lfloor \frac k2\rfloor} {k-s\choose s}(-\frac13)^s h_{k-s}=\frac1{3^k}-\frac 1{k 3^k}\)

\(\sum_{s=1}^{\lfloor \frac k2\rfloor} {k-s\choose s}(-1)^{s-1}3^{k-s} h_{k-s}=\frac{k-1}{k}\)
由此可以计算出$h_k$,比较有意思,在k不超过48时,$h_k$的绝对值都不大,但是超过48以后,会快速增加,比如50以内的近似值如下:

但是51项到60项如下:

ShuXueZhenMiHu 发表于 2022-6-25 02:42:45


页: [1] 2
查看完整版本: 数列 a(1)=1, a(n+1)=a(n)[1-a(n)/3], 求证 5/2<100a(100)<3