Galerkin-Legendre 谱方法求解分数阶微分方程的数值实现
一、由于任何偏微分方程初边值问题通过半离散之后都会转化为两点边值问题，以下仅以两点边值问题开始简要讨论Galerkin-Legendre 谱方法的数值实现过程
考虑如下分数阶两点边值问题
\begin{align}
&-(-\Delta)^{\frac{\alpha}{2}}u(x)=f(x), \quad a< x< b,\quad 1<\alpha\leq 2, \\
&u(a)=0,\quad u(b)=0,
\end{align}
其中
\begin{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}{2\cos(\frac{\alpha \pi}{2})},
\end{equation}
\begin{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)d\xi}{(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)d\xi}{(\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正交多项式，由以下三项递推公式确定
\begin{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}
\begin{align}
L_{0}(\hat{x})=1,\quad L_{1}(\hat{x})=\hat{x}.
\end{align}
选取空间$\mathbb{V}_N$中的基，对真解作如下插值
\begin{align}
u_N(x)=\sum\limits^{N-2}_{k=0}c_k\phi_k(x),
\end{align}
其中的$c_k$即为我们最终需要计算的量（注意：这里的插值并不是拉格朗日型的插值）。
我们可以把方程(1)的变分形式写成如下形式，finding $u_N\in\mathbb{V}_N$, 使得
\begin{align}
&\big( -(-\Delta)^{\frac{\alpha}{2}}u_N,v)=(f,v), \quad v\in\mathbb{V}_N,
\end{align}
其中, $(u,v)=\int^b_a u\bar{v}dx$,
\begin{align}
&\big( -(-\Delta)^{\frac{\alpha}{2}}u_N,v)=c_\alpha\Big( {}_{a}^{RL}\!D^{\alpha}_{x}u_N,v \Big) + c_\alpha\Big( {}_{x}^{RL}\!D^{\alpha}_{b}u_N,v \Big)
\end{align}
取$v=\phi_i(x),~~i=0,1,\cdots,N-2,$
\begin{align}
\Big( {}_{a}^{RL}\!D^{\alpha}_{x}u_N,\phi_i \Big)&=\sum\limits^{N-2}_{k=0}\Big( {}_{a}^{RL}\!D^{\alpha/2}_{x}\phi_k,{}_{x}^{RL}\!D^{\alpha/2}_{b}\phi_i \Big)c_k=\sum\limits^{N-2}_{k=0} (S_x)_{k,i}c_k,
\end{align}
\begin{align}
\Big( {}_{x}^{RL}\!D^{\alpha}_{b}u_N,\phi_i \Big)=\sum\limits^{N-2}_{k=0}\Big( {}_{x}^{RL}\!D^{\alpha/2}_{b}\phi_k,{}_{a}^{RL}\!D^{\alpha/2}_{x}\phi_i \Big)c_k=\sum\limits^{N-2}_{k=0} (S_y)_{k,i}c_k,
\end{align}
不难发现，$S_y=(S_x)^T$.
\begin{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}
其中
\begin{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)\nonumber\\ &= (\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\nonumber\\&= (\frac{b-a}{2})^{-\alpha} \frac{\Gamma(k+1)\Gamma(i+1)}{\Gamma(k-\alpha/2+1)\Gamma(i-\alpha/2+1)}\nonumber \\&\times\int^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\nonumber\\&= (\frac{b-a}{2})^{1-\alpha} \frac{\Gamma(k+1)\Gamma(i+1)}{\Gamma(k-\alpha/2+1)\Gamma(i-\alpha/2+1)}\nonumber \\&\times\int^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})d\hat{x}\nonumber\\&\thickapprox (\frac{b-a}{2})^{1-\alpha} \frac{\Gamma(k+1)\Gamma(i+1)}{\Gamma(k-\alpha/2+1)\Gamma(i-\alpha/2+1)}\nonumber \\&\times\sum\limits^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正交多项式
\begin{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)d\hat{x}\nonumber\\&\approx \frac{b-a}{2}\sum\limits^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正交多项式
格式写成矩阵的形式
\begin{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]内的表达形式即为
\begin{align}
u_N(x)=\sum\limits^{N-2}_{k=0}c_k\phi_k(x).
\end{align}
剩下的就是代码实现了...
