一、由于任何偏微分方程初边值问题通过半离散之后都会转化为两点边值问题，以下仅以两点边值问题开始简要讨论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}
剩下的就是代码实现了...