阿多米安分解法解半无限区域微分方程一例——气体在多孔介质中的扩散

参考文献: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 代码:

1
2
3
4
5
6
7
8
9
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 近似式:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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 近似式。

标签: 阿多米安分解法

添加新评论

captcha
请输入验证码