优化变分迭代法解Emden-Lane边值问题一例

参考文献：Eur. Phys. J. Plus (2017)132: 251

在前一篇文章中，给出了Emden-Lane边值问题的变分迭代法。

Emden-Lane边值问题：

\begin{equation} \begin{split} & u''(x)+\frac{\alpha}{x}u'(x)=f(x,u) \\ & u'(0)=0, au(1)+bu'(1)=c \end{split} \label{Emden-Lane} \end{equation}

本文计算一下$\alpha=1$的Emden-Lane边值问题。

由前一篇文章知，迭代格式：

\begin{equation} \begin{split} u_{n+1}(x)=&u_n(x)+\int_0^xs\ln\left(\frac{s}{x}\right)\left( u_n''(s)+\frac{1}{s}u_n'(s)-f(s,u_n(s)) \right) \mathrm ds \\ =&u_n(x)+\int_0^xs\ln\left(\frac{s}{x}\right)N(u_n(s)) \mathrm ds \end{split} \label{1iterpre} \end{equation}

在上式添加$u''(x)+\frac{1}{x}u'(x)=0$的解$c_1\ln x + c_2$以及辅助参数$h$，于是，迭代格式变为：

\begin{equation} u_{n+1}(x)=u_n(x)+h\int_0^xs\ln\left(\frac{s}{x}\right)N(u_n(s)) \mathrm ds+c_1\ln x + c_2 \label{1iter} \end{equation}

各$u_n$满足边界条件$u'(0)=0, au(1)+bu'(1)=c $，又由\eqref{1iter}式知，$u'_{n+1}(x)=u'_n(x)-h\int_0^x s N(u_n(s))/x \mathrm ds+c_1/x$，得

\begin{equation} c_1=0, c_2=h\int_0^1 s\left(\frac{b}{a}-\ln s\right)N(u_n(s)) \mathrm ds \label{c1c2} \end{equation}

代入\eqref{1iter}式得

\begin{equation} \begin{split} u_{n+1}(x)=&u_n(x)+h\int_0^xs\ln\left(\frac{s}{x}\right)N(u_n(s)) \mathrm ds+h\int_0^1 s\left(\frac{b}{a}-\ln s\right)N(u_n(s)) \mathrm ds\\ =& u_n(x)+h\int_0^x \left( \frac{b}{a}-\ln x \right) s N(u_n(s)) \mathrm ds+h\int_x^1 \left(\frac{b}{a}-\ln s\right) sN(u_n(s)) \mathrm ds \\ =& u_n(x)+h\int_0^1 G(x,s)[(su'_n)'-sf(s,u_n)]\mathrm ds \end{split} \label{1itercomput} \end{equation}

其中

\begin{equation} G(x,s)= \begin{cases} & \frac{b}{a}-\ln x , 0\le s \le x\\ & \frac{b}{a}-\ln s, x \le s \le 1 \end{cases} \label{G} \end{equation}

令

\begin{equation} u_n(x)=\sum_{k=0}^ny_n(x) \label{uy} \end{equation}

\eqref{1itercomput}式中的非线性项$f(s,u_n)$可展开称阿多米安多项式：

\begin{equation} f(s,u_n)=f(s,\sum_{k=0}^ny_n)=\sum_{k=0}^nA_k \label{fdom} \end{equation}

将\eqref{uy}式代入\eqref{1itercomput}式

\begin{equation} \sum_{k=0}^{n+1}y_k(x)=\sum_{k=0}^ny_k(x)+h\int_0^1 G(x,s)\left [\left(s\sum_{k=0}^ny'_k(s)\right)'-s\sum_{k=0}^nA_k\right ]\mathrm ds \label{yinterpre} \end{equation}

迭代格式为

\begin{equation} y_{n+1}(x)=h\int_0^1 G(x,s)\left [\left(s\sum_{k=0}^ny'_k(s)\right)'-s\sum_{k=0}^nA_k\right ]\mathrm ds \label{yinter} \end{equation}

辅助参数$h$由残差函数得到。

残差函数为：

\begin{equation} R_n(x,h)=u_n''(x,h)+\frac{1}{x}u_n'(x,h)-f(x,u_n(x,h)) \label{Rerror} \end{equation}

残差函数的平方和为：

\begin{equation} I_n(h)=\int_0^1 [R_n(x,h)]^2\mathrm dx \label{SRerror} \end{equation}

辅助参数$h$由下式得到：

\begin{equation} \frac{\partial I_n(h)}{\partial h}=\int_0^1 \frac{\partial [R_n(x,h)]^2}{\partial h}\mathrm dx=0 \label{h} \end{equation}

下面举个例子：

\begin{equation} \begin{split} & u''(x)+\frac{\alpha}{x}u'(x)+e^{u(x)}=0 \\ & u'(0)=0, u(1)=0 \end{split} \label{Exam1} \end{equation}

首先计算$f(u_n)=e^{u_n}$的阿多米安多项式，Mathematica代码为：

Do[A[k] = 
  Coefficient[
   Exp[y[0]]*Normal[Series[Exp[t], {t, 0, k}]] /. 
    t -> (Sum[y[j]*Exp[j*I*x], {j, 1, k}]), Exp[I*x], k], {k, 4}]
A[1]
A[2]
A[3]
A[4]

根据迭代格式：

\begin{equation} \begin{split} y_n(x)=&h\int_0^1 G(x,s)\left [\left(s\sum_{k=0}^ny'_k(s)\right)'+s\sum_{k=0}^nA_k(s)\right ]\mathrm ds \\ =& -h\ln x \int_0^x \left [\left(s\sum_{k=0}^ny'_k(s)\right)'+s\sum_{k=0}^nA_k(s)\right ]\mathrm ds \\ & -h\int_0^1 \ln s\left [\left(s\sum_{k=0}^ny'_k(s)\right)'+s\sum_{k=0}^nA_k(s)\right ]\mathrm ds \end{split} \label{yinterex1} \end{equation}

其中，

\begin{equation} G(x,s)= \begin{cases} & -\ln x , 0\le s \le x\\ & -\ln s, x \le s \le 1 \end{cases} \label{Gex1} \end{equation}

方程的解截断至$u_3$，Mathematica 代码为

A[0] = E^y[0]
y[0] = 0
y[1] = -h*Log[x]*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(x\)]\(\((D[
        x*D[y[0], x], x] + x*A[0])\) \[DifferentialD]x\)\) - h*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(x\), \(1\)]\(Log[
      x]*\((D[x*D[y[0], x], x] + x*A[0])\) \[DifferentialD]x\)\)
y[1] = Collect[Normal @ y[1], x]      
y[2] = -h*Log[x]*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(x\)]\(\((D[
        x*D[y[0] + y[1], x], x] + 
       x*\((A[0] + A[1])\))\) \[DifferentialD]x\)\) - h*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(x\), \(1\)]\(Log[
      x]*\((D[x*D[y[0] + y[1], x], x] + 
       x*\((A[0] + A[1])\))\) \[DifferentialD]x\)\)
y[2] = Collect[Normal @ y[2], x]
y[3] = -h*Log[x]*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(x\)]\(\((D[
        x*D[y[0] + y[1] + y[2], x], x] + 
       x*\((A[0] + A[1] + A[2])\))\) \[DifferentialD]x\)\) - h*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(x\), \(1\)]\(Log[
      x]*\((D[x*D[y[0] + y[1] + y[2], x], x] + 
       x*\((A[0] + A[1] + A[2])\))\) \[DifferentialD]x\)\)
y[3] = Collect[Normal @ y[3], x]
u[3] = Collect[y[0] + y[1] + y[2] + y[3], x, Simplify]

结果为：

\begin{equation} u_3(x,h)=-\frac{1}{768} h^3 x^6-\frac{3}{128} (h-2) h^2 x^4-\frac{1}{256} h \left(37 h^2-144 h+192\right) x^2+\frac{1}{384} h \left(65 h^2-234 h+288\right) \label{u3} \end{equation}

下面确定$h$。

直接由\eqref{Rerror}、\eqref{SRerror}、\eqref{h}三式无法计算出$h$，因为积分没有解析表达式。

将残差函数$R_3(x,h)$写成离散格式：

\begin{equation} R_n(i/N,h)=u_n''(i/N,h)+\frac{1}{x}u_n'(i/N,h)-f(x,u_n(i/N,h)) \label{Rerrordis} \end{equation}

残差函数的平方和为：

\begin{equation} I_n(h)=\sum_{i=0}^N [R_n(i/N,h)]^2 \label{SRerrorsum} \end{equation}

$h$由下式计算得到：

\begin{equation} \frac{\partial I_n(h)}{\partial h}=0 \label{hsol} \end{equation}

将区间$[0,1]$化出10等分，Mathematica计算代码：

R[3] = D[u[3], {x, 2}] + D[u[3], x]/x + Exp[u[3]] /. x -> i/10;
Ih[3] = Limit[D[u[3], {x, 2}] + D[u[3], x]/x + Exp[u[3]], x -> 0]^2 + 
   Sum[R[3]^2, {i, 1, 10}];
FindRoot[D[Ih[3], h] == 0, {h, 1}]

计算得$h=1.22414$。