差分方程求解问题

mathe

模仿7#代码有

1. ? 2/y+1/2*X+(b00+b01*X)*y+(b10+b11*X+b12*X^2)*y^2+(b20+b21*X+b22*X^2+b23*X^3)*y^3+O(y^4)
   
2. %13 = 2*y^-1 + 1/2*X + (b01*X + b00)*y + (b12*X^2 + b11*X + b10)*y^2 + (b23*X^3 + b22*X^2 + b21*X + b20)*y^3 + O(y^4)
   
3. ? subst(%,y,y/(1+y))
   
4. %14 = 2*y^-1 + (1/2*X + 2) + (b01*X + b00)*y + (b12*X^2 + (-b01 + b11)*X + (-b00 + b10))*y^2 + (b23*X^3 + (-2*b12 + b22)*X^2 + (b01 + (-2*b11 + b21))*X + (b00 + (-2*b10 + b20)))*y^3 + O(y^4)
   
5. ? subst(%,X,X+log(1+y))
   
6. %15 = 2*y^-1 + (1/2*X + 2) + (b01*X + (b00 + 1/2))*y + (b12*X^2 + (-b01 + b11)*X + (-b00 + (b01 + (b10 - 1/4))))*y^2 + (b23*X^3 + (-2*b12 + b22)*X^2 + (b01 + (-2*b11 + (2*b12 + b21)))*X + (b00 + (-3/2*b01 + (-2*b10 + (b11 + (b20 + 1/6))))))*y^3 + O(y^4)
   
7. ? %15-%13-2-1/%13
   
8. %16 = ((-b01 + 1/8)*X + (-b00 + (b01 - 1/4)))*y^2 + ((-2*b12 - 1/32)*X^2 + (5/4*b01 + (-2*b11 + 2*b12))*X + (5/4*b00 + (-3/2*b01 + (-2*b10 + (b11 + 1/6)))))*y^3 + O(y^4)
   
9. ? subst(%,b01,1/8)
   
10. %17 = (-b00 - 1/8)*y^2 + ((-2*b12 - 1/32)*X^2 + (-2*b11 + (2*b12 + 5/32))*X + (5/4*b00 + (-2*b10 + (b11 - 1/48))))*y^3 + O(y^4)
    
11. ? subst(%,b00,-1/8)
    
12. %18 = ((-2*b12 - 1/32)*X^2 + (-2*b11 + (2*b12 + 5/32))*X + (-2*b10 + (b11 - 17/96)))*y^3 + O(y^4)
    
13. ? subst(%,b12,-1/64)
    
14. %19 = ((-2*b11 + 1/8)*X + (-2*b10 + (b11 - 17/96)))*y^3 + O(y^4)
    
15. ? subst(%,b11,1/16)
    
16. %20 = (-2*b10 - 11/96)*y^3 + O(y^4)
    
17. ? subst(%,b10,-11/192)
    
18. %21 = O(y^4)
    

复制代码

所以得到结果为
2/y+1/2*X+(-1/8+1/8*X)*y+(-11/192+1/16*X-1/64*X^2)*y^2+...
这个结果和链接中是一致的

数学星空

根据楼上46#,47# mathe的分析结果我们试着推算出极限常数\(c_0\)渐近表达式：

\(m'=J(m)=d_{11}m+d_0+\frac{d_1}{m}+\frac{d_2}{m^2}+\frac{d_3}{m^3}+\frac{d_4}{m^4}+\frac{d_5}{m^5}+\frac{d_6}{m^6}+\frac{d_7}{m^7}\)

\(m(x)=e_{11}x+\frac{e_1}{x}+\frac{e_2}{x^2}+\frac{e_3}{x^3}+\frac{e_4}{x^4}+\frac{e_5}{x^5}+\frac{e_6}{x^6}+\frac{e_7}{x^7}\)

引用楼上的记号及结果

\(a_{n+1}=a_n+b_0+\frac{b_1}{a_n}+\frac{b_2}{a_n^2}+\frac{b_3}{a_n^3}+\frac{b_4}{a_n^4}+\frac{b_5}{a_n^5}+\frac{b_6}{a_n^6}+\frac{b_7}{a_n^7}+\frac{b_8}{a_n^8}\)

我们根据关系式

\(m(a_{n+1})=J(m(a_n))\)

渐近展开得到下面结果：

e11 = d0/b0,

e1 = (b0^2*b1 + 2*b0*b2 - 2*b1^2)/(2*b0^3),

e2 = (b0^4*b1 + 6*b0^3*b2 + 3*b0^2*b1^2 + 6*b0^2*b3 - 6*b1^3)/(12*b0^4),

e3 = (6*b0^5*b2 + 4*b0^4*b1^2 + 18*b0^4*b3 + 15*b0^2*b1^3 + 12*b0^3*b4 - 6*b0^2*b1*b3 - 12*b0^2*b2^2 + 36*b0*b1^2*b2 - 30*b1^4)/(36*b0^5),

e4 = -(6*b0^8*b1 - 30*b0^6*b1^2 - 180*b0^6*b3 - 180*b0^5*b1*b2 + 25*b0^4*b1^3 - 360*b0^5*b4 + 180*b0^4*b2^2 - 570*b0^2*b1^4 - 180*b0^4*b5 + 120*b0^3*b1*b4 + 360*b0^3*b2*b3 - 420*b0^2*b1^2*b3 + 60*b0^2*b1*b2^2 - 720*b0*b1^3*b2 + 780*b1^5)/(720*b0^6),

e5 = -(120*b0^9*b2 + 86*b0^8*b1^2 - 600*b0^7*b1*b2 - 30*b0^6*b1^3 - 1200*b0^7*b4 - 1260*b0^6*b1*b3 - 240*b0^6*b2^2 + 600*b0^5*b1^2*b2 + 215*b0^4*b1^4 - 1800*b0^6*b5 + 1800*b0^5*b2*b3 - 900*b0^4*b1^2*b3 + 900*b0^4*b1*b2^2 + 900*b0^3*b1^3*b2 - 5010*b0^2*b1^5 - 720*b0^5*b6 + 540*b0^4*b1*b5 + 1440*b0^4*b2*b4 + 720*b0^4*b3^2 - 1560*b0^3*b1^2*b4 - 720*b0^3*b1*b2*b3 - 720*b0^3*b2^3 - 1020*b0^2*b1^3*b3 + 4620*b0^2*b1^2*b2^2 - 9000*b0*b1^4*b2 + 6420*b1^6)/(3600*b0^7),

e6 = (600*b0^12*b1 - 3780*b0^10*b1^2 - 12600*b0^10*b3 - 16800*b0^9*b1*b2 - 1127*b0^8*b1^3 + 37800*b0^8*b1*b3 + 25200*b0^8*b2^2 + 12600*b0^7*b1^2*b2 - 5040*b0^6*b1^4 + 63000*b0^8*b5 + 67200*b0^7*b1*b4 + 25200*b0^7*b2*b3 - 5880*b0^6*b1^2*b3 - 50820*b0^6*b1*b2^2 + 14700*b0^5*b1^3*b2 - 62055*b0^4*b1^5 + 75600*b0^7*b6 - 75600*b0^6*b2*b4 - 37800*b0^6*b3^2 + 50400*b0^5*b1^2*b4 - 37800*b0^5*b1*b2*b3 + 25200*b0^5*b2^3 + 31500*b0^4*b1^3*b3 - 50400*b0^4*b1^2*b2^2 - 214200*b0^3*b1^4*b2 + 418320*b0^2*b1^6 + 25200*b0^6*b7 - 20160*b0^5*b1*b6 - 50400*b0^5*b2*b5 - 50400*b0^5*b3*b4 + 52920*b0^4*b1^2*b5 + 40320*b0^4*b1*b2*b4 + 26460*b0^4*b1*b3^2 + 75600*b0^4*b2^2*b3 + 15120*b0^3*b1^3*b4 - 246960*b0^3*b1^2*b2*b3 + 30240*b0^3*b1*b2^3 + 181440*b0^2*b1^4*b3 - 294840*b0^2*b1^3*b2^2 + 655200*b0*b1^5*b2 - 439740*b1^7)/(151200*b0^8),

e7 = (25200*b0^13*b2 + 22440*b0^12*b1^2 - 158760*b0^11*b1*b2 - 72870*b0^10*b1^3 - 176400*b0^11*b4 - 272160*b0^10*b1*b3 - 63000*b0^10*b2^2 + 2940*b0^9*b1^2*b2 - 22519*b0^8*b1^4 + 352800*b0^9*b1*b4 + 529200*b0^9*b2*b3 + 144900*b0^8*b1^2*b3 + 163800*b0^8*b1*b2^2 - 14700*b0^7*b1^3*b2 - 54705*b0^6*b1^5 + 529200*b0^9*b6 + 567000*b0^8*b1*b5 + 201600*b0^8*b2*b4 + 113400*b0^8*b3^2 + 63000*b0^7*b1^2*b4 - 718200*b0^7*b1*b2*b3 - 50400*b0^7*b2^3 + 280140*b0^6*b1^3*b3 + 322560*b0^6*b1^2*b2^2 - 345450*b0^5*b1^4*b2 - 1164135*b0^4*b1^6 + 529200*b0^8*b7 - 529200*b0^7*b2*b5 - 529200*b0^7*b3*b4 + 396900*b0^6*b1^2*b5 - 176400*b0^6*b1*b2*b4 + 529200*b0^6*b2^2*b3 + 352800*b0^5*b1^3*b4 - 970200*b0^5*b1^2*b2*b3 + 617400*b0^5*b1*b2^3 - 749700*b0^4*b1^3*b2^2 - 3748500*b0^3*b1^5*b2 + 5649840*b0^2*b1^7 + 151200*b0^7*b8 - 126000*b0^6*b1*b7 - 302400*b0^6*b2*b6 - 302400*b0^6*b3*b5 - 151200*b0^6*b4^2 + 312480*b0^5*b1^2*b6 + 289800*b0^5*b1*b2*b5 + 378000*b0^5*b1*b3*b4 + 453600*b0^5*b2^2*b4 + 453600*b0^5*b2*b3^2 + 30240*b0^4*b1^3*b5 - 1355760*b0^4*b1^2*b2*b4 - 633780*b0^4*b1^2*b3^2 - 50400*b0^4*b1*b2^2*b3 - 151200*b0^4*b2^4 + 861840*b0^3*b1^4*b4 - 1438920*b0^3*b1^3*b2*b3 + 1811880*b0^3*b1^2*b2^3 + 1862280*b0^2*b1^5*b3 - 6184080*b0^2*b1^4*b2^2 + 9437400*b0*b1^6*b2 - 5346180*b1^8)/(1058400*b0^9)]



d11 = 1,

d0 = 1,

d1 = b1/b0^2,

d2 = -b1*(b0^2 - 2*b1)/(2*b0^4),

d3 = b1*(2*b0^4 - 9*b0^2*b1 + 6*b1^2)/(6*b0^6),

d4 = -b1*(3*b0^6 - 22*b0^4*b1 + 36*b0^2*b1^2 - 12*b1^3)/(12*b0^8),

d5 = b1*(12*b0^8 - 125*b0^6*b1 + 350*b0^4*b1^2 - 300*b0^2*b1^3 + 60*b1^4)/(60*b0^10),

d6 = -b1*(20*b0^10 - 274*b0^8*b1 + 1125*b0^6*b1^2 - 1700*b0^4*b1^3 + 900*b0^2*b1^4 - 120*b1^5)/(120*b0^12),

d7 = b1*(120*b0^12 - 2058*b0^10*b1 + 11368*b0^8*b1^2 - 25725*b0^6*b1^3 + 24500*b0^4*b1^4 - 8820*b0^2*b1^5 + 840*b1^6)/(840*b0^14),


例1：对于\(a_{n+1}=a_n+3+\frac{3}{a_n}+\frac{1}{a_n^2}, a_1=1\)

\( b_0=3,b_1=3,b_2=1,b_3=0,b_4=0,b_5=0,b_6=0,b_7=0,b_8=0\) 代入上面计算结果得到

[e11 = d0/3, e1 = 5/18, e2 = 1/2, e3 = 239/324, e4 = 5227/6480, e5 = 7297/32400, e6 = -382271/226800, e7 = -12524549/3175200]

[d11 = 1, d0 = 1, d1 = 1/3, d2 = -1/18, d3 = -1/54, d4 = 7/324, d5 = -31/4860, d6 = -107/29160, d7 = 2833/612360, d8 = -641/459270]

分别计算

{n = 30, a[n] = 94.577173180254437058800379657718614764727625998154},
{n = 1000, a[n] = 3008.0454122322256151375839578413379920499777400494},
{n = 10000, a[n] = 30010.345913254673211205141938898372100854537725396},
{n = 20000, a[n] = 60011.038913478301168969286681242141053322422361195},
{n = 50000, a[n] = 150011.95510827412924639839086962443643695059717290}

分别代入：

m(x)=x/3 + 5/(18*x) + 1/(2*x^2) + 239/(324*x^3) + 5227/(6480*x^4) + 7297/(32400*x^5) - 382271/(226800*x^6) - 12524549/(3175200*x^7)

\(m_1=m(a_n)\)

得到：
{m1 = 31.528718222470931961040721399184452825742285880303, n = 30.},
{m1 = 1002.6818964776359736742726114064011875468126660498, n = 1000.},
{m1 = 10003.448647008180113922877077419762755563810842065, n = 10000.},
{m1 = 20003.679642455017249161974097023827907564665187251, n = 20000.},
{m1 = 50003.985037943102904049982074683246743056443826019, n = 50000.}

最后代入

\(\frac{b_1\ln(m_1)}{b_0}+c_0-b_0(m_1-n)=\ln(m_1)+c_0-3m_1+3n=0\)

算出极限常数\(c_0\)

\(n=30，      c_0=1.1352558473155026193075691268217083184824348896867\)
\(n=1000,    c_0=1.1352558473155037141953943468178728050585156348720\)
\(n=10000,  c_0=1.1352558473155037141953943477479294305468236191642\)
\(n=20000,  c_0=1.1352558473155037141953943477479294399504376590288\)
\(n=50000,  c_0=1.1352558473155037141953943477479294399873307675098\)


例2：对于\(a_{n+1}=a_n+2+\frac{1}{a_n},a_1=1\)

\( b_0=2,b_1=1,b_2=0,b_3=0,b_4=0,b_5=0,b_6=0,b_7=0,b_8=0\) 代入上面计算结果得到

[e11 = d0/2, e1 = 1/8, e2 = 11/96, e3 = 47/576, e4 = 371/11520, e5 = -2479/115200, e6 = -445883/9676800, e7 = 1664099/135475200],

[d11 = 1, d0 = 1, d1 = 1/4, d2 = -1/16, d3 = 1/192, d4 = 7/768, d5 = -39/5120, d6 = 173/61440, d7 = 229/573440, d8 = -5465/4128768]

分别计算

{n = 30, a[n] = 61.431432209680669889769949619696411575918064864861},
{n = 1000, a[n] = 2003.1776889787595201855576745536224560801916101574},
{n = 10000, a[n] = 20004.328408262068666673439189413258083518373043544},
{n = 20000, a[n] = 40004.674938335385378572226021542333781590837891074},
{n = 50000, a[n] = 100005.13305468207519580781135024383354791985937962}

分别代入：

m(x)=x/2 + 1/(8*x) + 11/(96*x^2) + 47/(576*x^3) + 371/(11520*x^4) - 2479/(115200*x^5) - 445883/(9676800*x^6) + 1664099/(135475200*x^7)

\(m_1=m(a_n)\)

得到：
{m1 = 30.717781610688269088439863150144092383876357179841, n = 30.},
{m1 = 1001.5889069187996813553058771447707422654784615311, n = 1000.},
{m1 = 10002.164210379968343003648581381986108478537734838, n = 10000.},
{m1 = 20002.337472292399101530415990149790165733522183032, n = 20000.},
{m1 = 50002.566528590984895252425829375081503605774792411, n = 50000.}

最后代入

\(\frac{b_1\ln(m_1)}{b_0}+c_0-b_0(m_1-n)=\frac{\ln(m_1)}{2}+c_0-2m_1+2n=0\)

算出极限常数\(c_0\)

\(n=30，      c_0=-0.27685762486257780836779529521519143374727922463262\)
\(n=1000,    c_0=-0.27685762486257653893643725183175348286104116846537\)
\(n=10000,  c_0=-0.27685762486257653893643725082357340651477821095325\)
\(n=20000,  c_0=-0.27685762486257653893643725082357339635784350297152\)
\(n=50000,  c_0=-0.27685762486257653893643725082357339631800579653841\)