Galerkin-Legendre 谱方法求解分数阶微分方程的数值实现

一、由于任何偏微分方程初边值问题通过半离散之后都会转化为两点边值问题，以下仅以两点边值问题开始简要讨论Galerkin-Legendre 谱方法的数值实现过程

考虑如下分数阶两点边值问题

egin{align}
&-(-Delta)^{frac{alpha}{2}}u(x)=f(x), quad  a< x< b,quad 1<alphaleq 2, \
&u(a)=0,quad u(b)=0,
end{align}

 其中

egin{equation}
-(-Delta)^{frac{alpha}{2}}u(x)(x)=c_alpha Big( {}_{a}^{RL}!D^{alpha}_{x}+{}_{x}^{RL}!D^{alpha}_{b} Big)u(x),qquad c_alpha=-frac{1}{2cos(frac{alpha pi}{2})},
end{equation}

egin{equation*}
{}_{a}^{RL}!D^{alpha}_{x}u(x)=frac{1}{Gamma(2-alpha)}frac{d^2}{d x^2}int_{a}^{x}frac{u(xi)dxi}{(x-xi)^{alpha-1}},quad
{}_{x}^{RL}!D^{alpha}_{b}u(x)=frac{1}{Gamma(2-alpha)}frac{d^2}{d x^2}int_{x}^{b}frac{u(xi)dxi}{(xi-x)^{alpha-1}},
end{equation*}

取基函数空间

$$mathbb{V}_N=span{  phi_0(x), phi_1(x),cdots, phi_{N-2}(x) },$$

且

$$phi_k(x)=L_k(hat{x})-L_{k+2}(hat{x}),~~hat{x}in[-1,1],quad x=frac{(b-a)hat{x}+(b+a)}{2}in[a,b],$$

其中

$L_k(hat{x})$为定义在[-1,1]内的Legendre正交多项式，由以下三项递推公式确定

egin{align}
& L_{k+1}(hat{x})=frac{2k+1}{k+1}hat{x}L_{k}(hat{x})-frac{k}{k+1}L_{k-1}(hat{x}), 
end{align}

egin{align}
 L_{0}(hat{x})=1,quad  L_{1}(hat{x})=hat{x}.
end{align}

选取空间$mathbb{V}_N$中的基，对真解作如下插值

egin{align}
u_N(x)=sumlimits^{N-2}_{k=0}c_kphi_k(x),
end{align}

其中的$c_k$即为我们最终需要计算的量（注意：这里的插值并不是拉格朗日型的插值）。

我们可以把方程(1)的变分形式写成如下形式，finding $u_Ninmathbb{V}_N$, 使得

egin{align}
&ig( -(-Delta)^{frac{alpha}{2}}u_N,v)=(f,v), quad  vinmathbb{V}_N,
end{align}

其中, $(u,v)=int^b_a uar{v}dx$, 

egin{align}
&ig( -(-Delta)^{frac{alpha}{2}}u_N,v)=c_alphaBig( {}_{a}^{RL}!D^{alpha}_{x}u_N,v  Big) +  c_alphaBig( {}_{x}^{RL}!D^{alpha}_{b}u_N,v  Big)
end{align}

取$v=phi_i(x),~~i=0,1,cdots,N-2,$

egin{align}
Big( {}_{a}^{RL}!D^{alpha}_{x}u_N,phi_i  Big)&=sumlimits^{N-2}_{k=0}Big( {}_{a}^{RL}!D^{alpha/2}_{x}phi_k,{}_{x}^{RL}!D^{alpha/2}_{b}phi_i   Big)c_k=sumlimits^{N-2}_{k=0} (S_x)_{k,i}c_k,
end{align}

egin{align}
Big( {}_{x}^{RL}!D^{alpha}_{b}u_N,phi_i  Big)=sumlimits^{N-2}_{k=0}Big( {}_{x}^{RL}!D^{alpha/2}_{b}phi_k,{}_{a}^{RL}!D^{alpha/2}_{x}phi_i   Big)c_k=sumlimits^{N-2}_{k=0} (S_y)_{k,i}c_k,
end{align}

不难发现，$S_y=(S_x)^T$.

egin{align}
Big( {}_{a}^{RL}!D^{alpha/2}_{x}phi_k,{}_{x}^{RL}!D^{alpha/2}_{b}phi_i   Big)&=Big( {}_{a}^{RL}!D^{alpha/2}_{x}L_k(hat{x}),{}_{x}^{RL}!D^{alpha/2}_{b}L_i(hat{x})   Big) -Big( {}_{a}^{RL}!D^{alpha/2}_{x}L_{k+2}(hat{x}),{}_{x}^{RL}!D^{alpha/2}_{b}L_i (hat{x})  Big)\
&~-Big( {}_{a}^{RL}!D^{alpha/2}_{x}L_{k}(hat{x}),{}_{x}^{RL}!D^{alpha/2}_{b}L_{i+2}  Big)+Big( {}_{a}^{RL}!D^{alpha/2}_{x}L_{k+2}(hat{x}),{}_{x}^{RL}!D^{alpha/2}_{b}L_{i+2} (hat{x})  Big),\&=D(k,i)-D(k+2,i)-D(k,i+2)
+D(k+2,i+2),\&quad k,i=0,1,cdots,N-2.end{align}

 其中

egin{align}   D(k,i)&=Big( {}_{a}^{RL}!D^{alpha/2}_{x}L_k(hat{x}),{}_{x}^{RL}!D^{alpha/2}_{b}L_i(hat{x})   Big) =(frac{b-a}{2})^{-alpha}Big( {}_{-1}^{RL}!D^{alpha/2}_{hat{x}}L_k(hat{x}),{}_{hat{x}}^{RL}!D^{alpha/2}_{1}L_i(hat{x})   Big) onumber\ &= (frac{b-a}{2})^{-alpha} int^b_a {}_{-1}^{RL}!D^{alpha/2}_{hat{x}}L_k(hat{x}){}_{hat{x}}^{RL}!D^{alpha/2}_{1}L_i(hat{x}) dx onumber\&= (frac{b-a}{2})^{-alpha} frac{Gamma(k+1)Gamma(i+1)}{Gamma(k-alpha/2+1)Gamma(i-alpha/2+1)} onumber \& imesint^b_a  (1+hat{x})^{-alpha/2} (1-hat{x})^{-alpha/2}  J_k^{alpha/2,-alpha/2}(hat{x})J_i^{-alpha/2,alpha/2}dx onumber\&= (frac{b-a}{2})^{1-alpha} frac{Gamma(k+1)Gamma(i+1)}{Gamma(k-alpha/2+1)Gamma(i-alpha/2+1)} onumber \& imesint^1_{-1}  (1+hat{x})^{-alpha/2} (1-hat{x})^{-alpha/2}  J_k^{alpha/2,-alpha/2}(hat{x})J_i^{-alpha/2,alpha/2}(hat{x})dhat{x} onumber\& hickapprox (frac{b-a}{2})^{1-alpha} frac{Gamma(k+1)Gamma(i+1)}{Gamma(k-alpha/2+1)Gamma(i-alpha/2+1)} onumber \& imessumlimits^M_{j=0}  w_j  J_k^{alpha/2,-alpha/2}(hat{x}_j)J_i^{-alpha/2,alpha/2}(hat{x}_j),end{align}

最后一步用到了Jacobi-Gauss-Lobatto数值积分，其中$J_i^{a,b}(hat{x})$表示定义在[-1,1]之间的Jacobi正交多项式，详细可以参考： Jacobi正交多项式

 egin{align}
F_i:=(f,phi_i(x))&=int^b_a f(x)phi_i(x)dx=frac{b-a}{2}int^1_{-1} f(x)phi_i(x)dhat{x} onumber\&approx frac{b-a}{2}sumlimits^M_{j=0} f(x_j) (L_i(hat{x})-L_{i+2}(hat{x}))hat{w}_j,quad i=0,1,cdots,N-2.
end{align}

这里用到了Legendre-Gauss-Lobatto数值积分，详细参考： Jacobi正交多项式

格式写成矩阵的形式

 egin{align}
&c_alpha( S_x+S_x^T )C=F,
end{align}

其中：

$$C=(c_0,c_1,cdots,c_{N-2})^T,qquad F=(F_0,F_1,cdots,F_{N-2})^T.$$

最后，我们的数值解在[a,b]内的表达形式即为

egin{align}
u_N(x)=sumlimits^{N-2}_{k=0}c_kphi_k(x).
end{align}

剩下的就是代码实现了...