参考文献:Applied Mathematics and Computation 118 (2001) 123-132
方程:
\begin{equation}
y''(x)+2x\frac{y'(x)}{\sqrt{1-\alpha y(x)}}=0,\quad 0<\alpha <1
\label{eq}
\end{equation}
边界条件:
\begin{equation}
y(0)=1,\quad \lim_{x\to \infty}y(x)=0
\label{bd}
\end{equation}
将\eqref{eq}积分,得
\begin{equation}
y(x)=1+Bx-2\int_0^x\int_0^xx\frac{y'(x)}{\sqrt{1-\alpha y(x)}}\mathrm dx\mathrm dx
\label{solform}
\end{equation}
设解为级数形式:$y(x)=\sum_{n=0}^{\infty}y_n(x)$,则\eqref{solform}变为
\begin{equation}
\sum_{n=0}^{\infty}y_n(x)=1+Bx-2\int_0^x\int_0^xx\sum_{n=0}^{\infty}A_n(x)\mathrm dx\mathrm dx
\label{Adosolform}
\end{equation}
其中$B=y'(0)$,待定。
级数$y(x)=\sum_{n=0}^{\infty}y_n(x)$中各项满足如下迭代关系:
\begin{equation}
\begin{split}
& y_0(x)=1 \\
& y_1(x)=Bx-2\int_0^x\int_0^xxA_0(x)\mathrm dx\mathrm dx \\
& y_{k+2}(x)=-2\int_0^x\int_0^xxA_{k+1}(x)\mathrm dx\mathrm dx,\quad k\ge 0
\end{split}
\label{Adoiter}
\end{equation}
下面求解阿多米安多项式。根据阿多米安多项式一种参数化算法,Mathematica 代码:
| | Do[A[k] = | | Coefficient[(Sum[yp[j]*Exp[j*I*x], {j, 0, k}])/Sqrt[ | | 1 - \[Alpha]*y[0]]*Normal[Series[1/Sqrt[1 - t], {t, 0, k}]] /. | | t -> ((\[Alpha]*(Sum[y[j]*Exp[j*I*x], {j, 1, k}]))/( | | 1 - \[Alpha]*y[0])), Exp[I*x], k], {k, 4}] | | Expand[A[1]] | | Expand[A[2]] | | Expand[A[3]] | | Expand[A[4]] |
|
与文献所得结果一致。
下面计算级数解各项:
\begin{equation}
\begin{split}
& y_0(x)=1 \\
& y_1(x)=Bx-2\int_0^x\int_0^xxA_{0}(x)\mathrm dx\mathrm dx =Bx \\
& y_{2}(x)=-2\int_0^x\int_0^xxA_{1}(x)\mathrm dx\mathrm dx=-\frac{B x^3}{3 \sqrt{1-\alpha }} \\
& y_{3}(x)=-2\int_0^x\int_0^xxA_{2}(x)\mathrm dx\mathrm dx=\frac{\alpha B^2 x^4}{12 \sqrt{1-\alpha } (\alpha -1)}-\frac{B x^5}{10 (\alpha -1)} \\
& y_{4}(x)=-2\int_0^x\int_0^xxA_{3}(x)\mathrm dx\mathrm dx=-\frac{3 \alpha ^2 B^3 x^5}{80 (1-\alpha )^{5/2}}+\frac{\alpha B^2 x^6}{15 (1-\alpha )^2}-\frac{B x^7}{42 (1-\alpha )^{3/2}} \\
& y_{5}(x)=-\frac{\alpha ^3 B^4 x^6}{48 (1-\alpha )^{7/2}}+\frac{7 \alpha ^2 B^3 x^7}{144 (1-\alpha )^3}-\frac{13 \alpha B^2 x^8}{420 (1-\alpha )^{5/2}}+\frac{B x^9}{216 (1-\alpha )^2}
\end{split}
\label{Adosol}
\end{equation}
Mathematica 代码:
| 1 | | 2 | | 3 | | 4 | | 5 | | 6 | | 7 | | 8 | | 9 | | 10 | | 11 | | 12 | | 13 | | 14 | | 15 | | 16 | | 17 | | 18 | | 19 | | 20 | | 21 | | 22 | | 23 | | 24 | | 25 | | 26 | | 27 | | 28 | | 29 | | 30 | | 31 | | 32 | | 33 | | 34 | | 35 | | 36 |
| | In[1]:= A[0] = yp[0]/Sqrt[1 - \[Alpha] y[0]] | | Do[A[k] = | | Coefficient[(Sum[yp[j]*Exp[j*I*x], {j, 0, k}])/Sqrt[ | | 1 - \[Alpha]*y[0]]*Normal[Series[1/Sqrt[1 - t], {t, 0, k}]] /. | | t -> ((\[Alpha]*(Sum[y[j]*Exp[j*I*x], {j, 1, k}]))/( | | 1 - \[Alpha]*y[0])), Exp[I*x], k], {k, 4}] | | <p>Out[1]= yp[</p> | | <p>Out[9]= 1</p> | | <p>Out[10]= 0</p> | | <p>Out[11]= B x<br /> | | In[13]:= yp[1] = D[y[1], x]</p> | | <p>Out[13]= B</p> | | <p>In[14]:= y[2] = -2<em>Integrate[Integrate[x</em>A[1], {x, 0, x}], {x, 0, x}]</p> | | <p>Out[14]= -((B x^3)/(3 Sqrt[1 - [Alpha]]))<br /> | | In[15]:= yp[2] = D[y[2], x]<br /> | | y[3] = -2<em>Integrate[Integrate[x</em>A[2], {x, 0, x}], {x, 0, x}]</p> | | <p>Out[15]= -((B x^2)/Sqrt[1 - [Alpha]])</p> | | <p>Out[16]= -2 (-((B x^5)/(20 (1 - [Alpha]))) + (B^2 x^4 [Alpha])/(<br /> | | 24 (1 - [Alpha])^(3/2)))</p> | | <p>In[17]:= Expand[-2 (-((B x^5)/(20 (1 - [Alpha]))) + (<br /> | | B^2 x^4 [Alpha])/(24 (1 - [Alpha])^(3/2)))]</p> | | <p>Out[17]= (B x^5)/(10 (1 - [Alpha])) - (B^2 x^4 [Alpha])/(<br /> | | 12 (1 - [Alpha])^(3/2))<br /> | | In[18]:= yp[3] = D[y[3], x]<br /> | | y[4] = -2<em>Integrate[Integrate[x</em>A[3], {x, 0, x}], {x, 0, x}]</p> | | <p>Out[18]= -2 (-((B x^4)/(4 (1 - [Alpha]))) + (B^2 x^3 [Alpha])/(<br /> | | 6 (1 - [Alpha])^(3/2)))</p> | | <p>Out[19]= -2 ((B x^7)/(84 (1 - [Alpha])^(3/2)) - (B^2 x^6 [Alpha])/(<br /> | | 30 (1 - [Alpha])^2) + (3 B^3 x^5 [Alpha]^2)/(<br /> | | 160 (1 - [Alpha])^(5/2)))</p> | | <p>In[20]:= Expand[-2 ((B x^7)/(84 (1 - [Alpha])^(3/2)) - (<br /> | | B^2 x^6 [Alpha])/(30 (1 - [Alpha])^2) + (3 B^3 x^5 [Alpha]^2)/(<br /> | | 160 (1 - [Alpha])^(5/2)))]</p> | | <p>Out[20]= -((B x^7)/(42 (1 - [Alpha])^(3/2))) + (B^2 x^6 [Alpha])/(<br /> | | 15 (1 - [Alpha])^2) - (3 B^3 x^5 [Alpha]^2)/(<br /> | | 80 (1 - [Alpha])^(5/2))</p> |
|
文章Applied Mathematics and Computation 118 (2001) 123-132将解只保留到$x$的六阶项:
\begin{equation}
y(x)=x^5 \left(\frac{B}{10 (1-\alpha )}-\frac{3 \alpha ^2 B^3}{80 (1-\alpha )^{5/2}}\right)-\frac{\alpha B^2 x^4}{12 (1-\alpha )^{3/2}}+x^6 \left(\frac{\alpha B^2}{15 (1-\alpha )^2}-\frac{\alpha ^3 B^4}{48 (1-\alpha )^{7/2}}\right)-\frac{B x^3}{3 \sqrt{1-\alpha }}+B x+1
\label{Adosol6}
\end{equation}
[2/2]Pade 近似式:
| | In[23]:= PadeApproximant[ys, {x, 0, {2, 2}}] | | <p>Out[23]= (1 + (B x (-4 + 5 [Alpha]))/(4 (-1 + [Alpha])) + (<br /> | | x^2 (-4 + 4 [Alpha] + 3 B^2 Sqrt[1 - [Alpha]] [Alpha]))/(<br /> | | 12 Sqrt[1 - [Alpha]] (-1 + [Alpha])))/(1 + x^2/(<br /> | | 3 Sqrt[1 - [Alpha]]) + (B x [Alpha])/(4 (-1 + [Alpha])))</p> | | <p>In[24]:= Simplify[(<br /> | | 1 + (B x (-4 + 5 [Alpha]))/(4 (-1 + [Alpha])) + (<br /> | | x^2 (-4 + 4 [Alpha] + 3 B^2 Sqrt[1 - [Alpha]] [Alpha]))/(<br /> | | 12 Sqrt[1 - [Alpha]] (-1 + [Alpha])))/(<br /> | | 1 + x^2/(3 Sqrt[1 - [Alpha]]) + (B x [Alpha])/(4 (-1 + [Alpha])))]</p> | | <p>Out[24]= (-12 (1 - [Alpha])^(3/2) +<br /> | | 3 B x Sqrt[1 - [Alpha]] (-4 + 5 [Alpha]) +<br /> | | x^2 (-4 + (4 +<br /> | | 3 B^2 Sqrt[1 - [Alpha]]) [Alpha]))/(-12 (1 - [Alpha])^(<br /> | | 3/2) + 4 x^2 (-1 + [Alpha]) + 3 B x Sqrt[1 - [Alpha]] [Alpha])</p> |
|
根据边界条件,$\lim_{x\to \infty}y(x)=0$,得Pade 近似式中分子最高阶项系数为0,即$\alpha \left(3 \sqrt{1-\alpha } B^2+4\right)-4=0$,得
\begin{equation}
B=-\frac{2(1-\alpha)^{1/4}}{\sqrt{3\alpha}}
\label{B}
\end{equation}
与文献结果一致。
同理可计算[3/3]Pade 近似式。