Caputo 分数阶一维问题基于 L2-1σ逼近的快速差分方法(附Matlab源代码)

Caputo 分数阶一维问题基于 L2-1σ逼近的快速差分方法

Caputo 分数阶一维问题基于 L2-1${}_{\sigma }$ 逼近的空间二阶方法(附Matlab程序)

考虑如下时间分数阶慢扩散方程初边值问题
$$\{\begin{array}{l}{}_0^C{D}_t^{\alpha }u(x,t)=u_{xx}(x,t)+f(x,t),x\in (0,L),t\in (0,T]\\ u(x,0)=\phi (x),x\in (0,L)\\ u(0,t)=\mu (t),u(L,t)=\nu (t),t\in [0,T]\end{array}$$
其中 $\alpha \in (0,1),f,\phi ,\mu ,\nu$ 为已知函数, 且 $\phi (0)=\mu (0),\phi (L)=\nu (0)$.

设解函数 $u\in C^{(4,2)}([0,L]\times [0,T])$.

本文给出求解上述时间分数阶慢扩散方程初边值问题基于 L2-1${}_{\sigma }$ 逼近的快速差分方法.

差分格式的建立

记
$$\sigma =1-\frac{\alpha }{2},t_{n-1+\sigma }=(n-1+\sigma )\tau ,f_i^{n-1+\sigma }=f\left(x_i,t_{n-1+\sigma }\right)$$

在点 $\left(x_i,t_{n-1+\sigma }\right)$ 处考虑微分方程, 得到
$${}_0^C{D}_t^{\alpha }u\left(x_i,t_{n-1+\sigma }\right)=u_{xx}\left(x_i,t_{n-1+\sigma }\right)+f_i^{n-1+\sigma },1⩽i⩽M-1,1⩽n⩽N.$$
对上式中时间分数阶导数应用参考文献中快速逼近的理论, 可得
$$\{\begin{array}{c}{}_0^C{D}_t^{\alpha }u\left(x_i,t_{n-1+\sigma }\right)=\frac{1}{\Gamma (1-\alpha )}{\sum }_{l=1}^{N_{\exp }}{\omega }_l{F}_{l,i}^n+d_0^{(1,\alpha )}\left(U_i^n-U_i^{n-1}\right)+O\left({\tau }^{3-\alpha }+ϵ\right),\\ 1⩽i⩽M-1,1⩽n⩽N,\\ F_{l,i}^1=0,1⩽l⩽N_{\exp },1⩽i⩽M-1,\\ F_{l,i}^n={\mathrm{e}}^{-s_l\tau }{F}_{l,i}^{n-1}+A_l\left(U_i^{n-1}-U_i^{n-2}\right)+B_l\left(U_i^n-U_i^{n-1}\right),\\ 1⩽l⩽N_{\exp },1⩽i⩽M-1,2⩽n⩽N.\end{array}$$
对空间二阶导数应用二阶中心差商离散, 可得
$$\begin{array}{rl}{u}_{xx}\left(x_i,t_{n-1+\sigma }\right)& =\sigma u_{xx}\left(x_i,t_n\right)+(1-\sigma )u_{xx}\left(x_i,t_{n-1}\right)+O\left({\tau }^2\right)\\ & =\sigma {\delta }_x^2{U}_i^n+(1-\sigma ){\delta }_x^2{U}_i^{n-1}+O\left(h^2\right)+O\left({\tau }^2\right)\end{array}$$
于是得到
$$\{\begin{array}{l}\frac{1}{\Gamma (1-\alpha )}{\sum }_{l=1}^{N_{\text{exp }}}{\omega }_l{F}_{l,i}^n+d_0^{(1,\alpha )}\left(U_i^n-U_i^{n-1}\right)=\sigma {\delta }_x^2{U}_i^n+(1-\sigma ){\delta }_x^2{U}_i^{n-1}+f_i^{n-1+\sigma }+{\left(r_7\right)}_i^n,\\ F_{l,i}^1=0,1⩽i⩽M-1,1⩽n⩽N,\\ F_{l,i}^n={\mathrm{e}}^{-s_l\tau }{F}_{l,i}^{n-1}+A_l\left(U_i^{n-1}-U_i^{n-2}\right)+B_l\left(U_i^n-U_i^{n-1}\right),\\ 1⩽l⩽N_{\exp },1⩽i⩽M-1,2⩽n⩽N,\end{array}$$
且存在正常数 $c_7$ 使得
$$\left|{\left(r_7\right)}_i^n\right|⩽c_7\left({\tau }^2+h^2+ϵ\right),1⩽i⩽M-1,1⩽n⩽N.$$
消去中间变量 $\left\{F_{l,i}^n\right\}$ 可得
$$\begin{array}{rl}\sum _{k=0}^{n-1}{d}_k^{(n,\alpha )}\left(U_i^{n-k}-U_i^{n-k-1}\right)& =\sigma {\delta }_x^2{U}_i^n+(1-\sigma ){\delta }_x^2{U}_i^{n-1}+f_i^{n-1+\sigma }+{\left(r_7\right)}_i^n\\ 1& ⩽i⩽M-1,1⩽n⩽N\end{array}$$
注意到初边值条件, 有
$$\{\begin{array}{ll}{U}_i^0=\phi \left(x_i\right),& 1⩽i⩽M-1,\\ U_0^n=\mu \left(t_n\right),& U_M^n=\nu \left(t_n\right),0⩽n⩽N\end{array}$$

在上式中略去小量项 ${\left(r_7\right)}_i^n$, 并用数值解 $u_i^n$ 代替精确解 $U_i^n$, 可得求解问题的如下差分格式:
$$\{\begin{array}{l}\frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^{N_{\exp }}{\omega }_l{F}_{l,i}^n+d_0^{(1,\alpha )}\left(u_i^n-u_i^{n-1}\right)=\sigma {\delta }_x^2{u}_i^n+(1-\sigma ){\delta }_x^2{u}_i^{n-1}+f_i^{n-1+\sigma },\\ F_{l,i}^1=0,1⩽l⩽N_{\exp },1⩽i⩽M-1,\\ F_{l,i}^n={\mathrm{e}}^{-s_l\tau }{F}_{l,i}^{n-1}+A_l\left(u_i^{n-1}-u_i^{n-2}\right)+B_l\left(u_i^n-u_i^{n-1}\right),\\ 1⩽l⩽N_{\exp },1⩽i⩽M-1,2⩽n⩽N,\\ u_i^0=\phi \left(x_i\right),1⩽i⩽M-1,\\ u_0^n=\mu \left(t_n\right),u_M^n=\nu \left(t_n\right),0⩽n⩽N.\end{array}$$
其中
$$d_0^{(1,\alpha )}=\frac{{\sigma }^{1-\alpha }{\tau }^{-\alpha }}{\Gamma (2-\alpha )}$$

$$A_l={\int }_0^1\left(\frac{3}{2}-\xi \right){\mathrm{e}}^{-s_l(\sigma +1-\xi )\tau }\mathrm{d}\xi ,B_l={\int }_0^1\left(\xi -\frac{1}{2}\right){\mathrm{e}}^{-s_l(\sigma +1-\xi )\tau }\mathrm{d}\xi$$

分数阶微分方程数值算例

数值算例

考虑问题:
$${}_0^C{D}_t^{\alpha }u(x,t)=\frac{{\partial }^{2}u}{\partial x^2}(x,t)+f(x,t),(x,t)\in (0,1)\times (0,1),$$

其中 $0<\alpha <1$, 右端源项 $f(x,t)$ 和相应的初边值条件由真解 $u(x,t)=\sin (x)t^2$ 确定.

源项和初边值条件

由真解 $u(x,t)=\sin (x)t^2.$
则有
$$\{\begin{array}{l}{u}_x=t^2\cos x\\ u_{xx}=-t^2\sin x\end{array}\{\begin{array}{l}{u}_t=2t\sin x\\ u_{tt}=2\sin x.\end{array}$$

将真解对应部分带入微分方程从而反解右端源项 $f(x,t)$

$$\begin{array}{rl}{}_0^c{D}_t^{\alpha }u(x,t)& =\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{2s\sin x}{(t-s{)}^{\alpha }}ds\\ & =\frac{-2\sin x}{\Gamma (1-\alpha )}{\int }_0^t(t-s-t)(t-s{)}^{-\alpha }ds\\ & =\frac{-2\sin x}{\Gamma (1-\alpha )}{\int }_0^t\left[(t-s{)}^{1-\alpha }-t(t-s{)}^{-\alpha }\right]ds\\ & =\frac{-2\sin x}{\Gamma (1-\alpha )}\left[{\int }_0^t-(t-s{)}^{1-\alpha }d(t-s)+{\int }_0^{1}t(t-s{)}^{-\alpha }d(t-s)\right]\\ & =\frac{-2\sin x}{\Gamma (1-\alpha )}\left[{\frac{1}{2-\alpha }(t-s{)}^{2-\alpha }|}_t^0+{\frac{t}{1-\alpha }(t-s{)}^{1-\alpha }|}_0^t\right]\\ & =\frac{-2\sin x}{\Gamma (1-\alpha )}\left[\frac{t^{2-\alpha }}{2-\alpha }+\frac{t}{1-\alpha }\cdot \left(-t^{1-\alpha }\right)\right]\\ & =\frac{-2\sin x}{\Gamma (1-\alpha )}\cdot \frac{-t^{2-\alpha }}{(2-\alpha )(1-\alpha )}\\ & =\frac{2\sin x}{\Gamma (3-\alpha )}{t}^{2-\alpha }\end{array}$$

由于$${}_0^C{D}_t^{\alpha }u(x,t)=\frac{{\partial }^{2}u}{\partial x^2}(x,t)+f(x,t)$$
于是$$\begin{array}{rl}f(x,t)& ={}_0^C{D}_t^{\alpha }u(x,t)-\frac{{\partial }^{2}u}{\partial x^2}(x,t)\\ & =\frac{2\sin x}{\Gamma (3-\alpha )}{t}^{2-\alpha }+t^2\sin x.\end{array}$$

综上所述, 完整的数值算例为:

$$\{\begin{array}{l}{}_0^C{D}_t^{\alpha }u(x,t)=u_{xx}(x,t)+\frac{2\sin x}{\Gamma (3-\alpha )}{t}^{2-\alpha }+t^2\sin x,x\in (0,L),t\in (0,T]\\ u(x,0)=0,x\in (0,L)\\ u(0,t)=0,u(L,t)=t^2\sin L,t\in [0,T]\end{array}$$

代数系统

由于 $$F_{l,i}^n={\mathrm{e}}^{-s_l\tau }{F}_{l,i}^{n-1}+A_l\left(u_i^{n-1}-u_i^{n-2}\right)+B_l\left(u_i^n-u_i^{n-1}\right)$$ 中含有 $u_i^n$, 故考虑将 $F_{l,i}^n$ 拆成
$$F_{l,i}^n={\hat{F}}_{l,i}^n+B_l\left(u_i^n-u_i^{n-1}\right)$$
其中 $${\hat{F}}_{l,i}^n={\mathrm{e}}^{-s_l\tau }{F}_{l,i}^{n-1}+A_l\left(u_i^{n-1}-u_i^{n-2}\right)$$
根据差分格式:
$$\frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^{N_{\exp }}{\omega }_l{F}_{l,i}^n+d_0^{(1,\alpha )}\left(u_i^n-u_i^{n-1}\right)=\sigma {\delta }_x^2{u}_i^n+(1-\sigma ){\delta }_x^2{u}_i^{n-1}+f_i^{n-1+\sigma }$$
于是有
$$\frac{1}{\Gamma (1-\alpha )}\sum _{i=1}^{N_{\text{exp }}}{w}_i{\hat{F}}_{l,i}^n+\frac{1}{\Gamma (1-\alpha )}\sum _{i=1}^{N_{exp}}{\omega }_i{B}_l\left(u_i^n-u_i^{n-1}\right)+d_0^{(1,\alpha )}\left(u_i^n-u_i^{n-1}\right)$$
$$=\sigma {\delta }_x^2{u}_i^n+(1-\sigma ){\delta }_x^2{u}_i^{n-1}+f_i^{n-1+\sigma }$$

得到如下代数系统:
$$\mathbf{A}\left(\begin{array}{c}{u}_1^n\\ u_2^n\\ ⋮\\ u_{M-2}^n\\ u_{M-1}^n\end{array}\right)=\mathbf{b}$$
其中
$$\mathbf{A}=\left(\begin{array}{cccc}\frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^{N_{exp}}{\omega }_l{B}_l+d_0^{(1,\alpha )}+\frac{2\sigma }{h^2}& -\frac{\sigma }{h^2}& & \\ -\frac{\sigma }{h^2}& \frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^{N_{exp}}{\omega }_l{B}_l+d_0^{(1,\alpha )}+\frac{2\sigma }{h^2}& -\frac{\sigma }{h^2}& \\ & \ddots & \ddots & \\ & & \frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^{N_{exp}}{\omega }_l{B}_l+d_0^{(1,\alpha )}+\frac{2\sigma }{h^2}& -\frac{\sigma }{h^2}\end{array}\right)$$

$$\begin{array}{rl}\mathbf{b}& =\left(\begin{array}{cccc}\frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^{N_{exp}}{\omega }_l{B}_l+d_0^{(1,\alpha )}-\frac{2(1-\sigma )}{h^2}& \frac{1-\sigma }{h^2}& & \\ \frac{1-\sigma }{h^2}& \frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^{N_{exp}}{\omega }_l{B}_l+d_0^{(1,\alpha )}-\frac{2(1-\sigma )}{h^2}& \frac{1-\sigma }{h^2}& \\ & \ddots & \ddots & \\ & & \frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^{N_{exp}}{\omega }_l{B}_l+d_0^{(1,\alpha )}-\frac{2(1-\sigma )}{h^2}& \frac{1-\sigma }{h^2}\end{array}\right)\left(\begin{array}{c}{u}_1^{n-1}\\ u_2^{n-1}\\ ⋮\\ u_{M-2}^{n-1}\\ u_{M-1}^{n-1}\end{array}\right)\\ & +\frac{1-\sigma }{h^2}\cdot \left(\begin{array}{c}{u}_0^{n-1}\\ 0\\ ⋮\\ 0\\ u_{M-1}^n\end{array}\right)+\frac{\sigma }{h^2}\cdot \left(\begin{array}{c}{u}_0^n\\ 0\\ ⋮\\ 0\\ u_M^n\end{array}\right)+\frac{{\hat{a}}_0^{(\alpha )}}{\Gamma (1-\alpha )}\cdot \left(\begin{array}{c}{u}_1^{n-1}\\ u_2^{n-1}\\ ⋮\\ u_{M-2}^{n-1}\\ u_{M-1}^{n-1}\end{array}\right)-\frac{1}{\Gamma (1-\alpha )}\sum _{l=1}^{N_{exp}}{\omega }_l{\hat{F}}_{l,i}^n+f_i^{n-1+\sigma }\end{array}$$

相关积分计算:
$$\begin{array}{rl}{A}_l& ={\int }_0^1\left(\frac{3}{2}-\xi \right){\mathrm{e}}^{-s_l(\sigma +1-\xi )\tau }\mathrm{d}\xi \\ & ={\int }_0^1\frac{3}{2}{e}^{-s_l(\sigma +1-\xi )\tau }d\xi -{\int }_0^1\xi e^{-s_l(\sigma +1-\xi )\tau }d\xi \\ & =\frac{3}{2s_l\tau }\left[e^{-s_l\tau \sigma }-e^{-s_l\tau (\sigma +1)}\right]-e^{-s_l\left(\sigma +1\right)\tau }\cdot \frac{1}{s_l\tau }\left[e^{s_l\tau }-\frac{1}{s_{\tau }\tau }{e}^{s_l\tau }+\frac{1}{s_l\tau }\right]\\ & =\frac{1}{s_l\tau }{e}^{-s_l\tau \sigma }\left[\frac{1}{2}+\frac{1}{s_l\tau }-\left(\frac{1}{s_l\tau }+\frac{3}{2}\right)e^{-s_l\tau }\right]\end{array}$$
$$\begin{array}{rl}{B}_l& ={\int }_0^1\left(\xi -\frac{1}{2}\right){\mathrm{e}}^{-s_l(\sigma +1-\xi )\tau }\mathrm{d}\xi \\ & ={\int }_0^1-\frac{1}{2}{e}^{-s_l(\sigma +1-\xi )\tau }d\xi +{\int }_0^1\xi e^{-s_l(\sigma +1-\xi )\tau }d\xi \\ & =-\frac{1}{2s_l\tau }\left[e^{-s_l\tau \sigma }-e^{-s_l\tau (\sigma +1)}\right]+e^{-s_l\left(\sigma +1\right)\tau }\cdot \frac{1}{s_l\tau }\left[e^{s_l\tau }-\frac{1}{s_{\tau }\tau }{e}^{s_l\tau }+\frac{1}{s_l\tau }\right]\\ & =\frac{1}{s_l\tau }{e}^{-s_l\tau \sigma }\left[\frac{1}{2}-\frac{1}{s_l\tau }+\left(\frac{1}{s_l\tau }+\frac{1}{2}\right)e^{-s_l\tau }\right]\end{array}$$

数值结果

时间方向参数选取:

时间方向数值结果:

空间方向参数选取:

空间方向数值结果:

程序运行时间:

时间步长为 0.000100 空间步长为 0.012500 时快速 L2-1${}_{\sigma }$ 插值逼近历时 4.232766 秒.

误差Error =6.6569e-06

而相同的网格剖分及精度下, 传统的基于 L2-1${}_{\sigma }$ 逼近需历时 166.610357 秒.

误差6.9756e-07

源代码

主程序:

clc,clear
tic
%% 初始化定解问题
alpha=0.2;
sigma=1-alpha/2;
epsilon=1e-8;
tau_hat=1e-6;
tau=1/10000;   T_a=0;      T_b=1;  
h=1/80;     L_a=0;      L_b=1;                
M=(L_b-L_a)/h;  
N=(T_b-T_a)/tau;       
x=L_a:h:L_b;        t=T_a:tau:T_b; 

u=zeros(M+1,N+1);          %   @u 初始化数值解
u(:,1)=f_ic(x,T_a,0,0);    %   定义初值条件
u(1,:)=f_bc(L_a,t,0,0);    %   定义左边界条件
u(M+1,:)=f_bc(L_b,t,0,0);  %   定义右边界条件

[xs,ws,nexp] = sumofexpappr2new(alpha,epsilon,tau_hat,T_b);
%% 两层格式（启动层）
n=2;
F_l=zeros(M-1,nexp);
% 创建代数系统
% 系数矩阵 Coefficient_Matrix_A
Coefficient_Matrix_A=toeplitz([d(0,1,alpha,sigma,tau)+2*sigma/h^2,-sigma/h^2,zeros(M-3,1)']);
% 右端 Right_Term_B
Right_Term_B=f_source(x(2:M),t(n-1)+sigma*tau,alpha)'...
    -1/gamma(1-alpha)*F_l*ws...
    +toeplitz([d(0,1,alpha,sigma,tau)-2*(1-sigma)/h^2,(1-sigma)/h^2,zeros(M-3,1)'])*u(2:M,n-1)...
    +(1-sigma)/h^2*[u(1,n-1);zeros(M-3,1);u(M+1,n-1)]...
    +sigma/h^2*[u(1,n);zeros(M-3,1);u(M+1,n)];
% 求解迭代系统（由0层求第1层的值）
u(2:M,n)=Coefficient_Matrix_A\Right_Term_B;
fprintf('进程:\t%d/%d\n',n,N+1)

%% 三层格式
% A_l=3./(2*xs*tau).*(exp(1).^(-xs*tau*sigma)-exp(1).^(-xs*tau*(sigma+1)))...
%     -exp(1).^(-xs*(sigma+1)*tau)./(xs*tau).*(exp(1.^(xs*tau)-1./(xs*tau).*exp(1).^(xs*tau)+1./(xs*tau)));
% B_l=-1./(2*xs*tau).*(exp(1).^(-xs*tau*sigma)-exp(1).^(-xs*tau*(sigma+1)))...
%     +exp(1).^(-xs*(sigma+1)*tau)./(xs*tau).*(exp(1.^(xs*tau)-1./(xs*tau).*exp(1).^(xs*tau)+1./(xs*tau)));
% 使用上述未化简的 A_l、B_l 会出现较大误差, 考虑是因为较大舍入误差造成
A_l=1./(xs*tau).*exp(1).^(-xs*tau*sigma).*(1/2+1./(xs*tau)-(1./(xs*tau)+3/2).*exp(1).^(-xs*tau));
B_l=1./(xs*tau).*exp(1).^(-xs*tau*sigma).*(1/2-1./(xs*tau)+(1./(xs*tau)+1/2).*exp(1).^(-xs*tau));
for n=3:N+1
    % 创建代数系统
    F_l=F_l*diag(exp(1).^(-xs*tau))+((u(2:M,n-1)-u(2:M,n-2)))*A_l';
    % 系数矩阵 Coefficient_Matrix_A
    Coefficient_Matrix_A=toeplitz([1/gamma(1-alpha)*B_l'*ws+d(0,1,alpha,sigma,tau)+2*sigma/h^2,-sigma/h^2,zeros(M-3,1)']);
    % 右端 Right_Term_B1
    Right_Term_B1=f_source(x(2:M),t(n-1)+sigma*tau,alpha)'...
        -1/gamma(1-alpha)*F_l*ws...
        +toeplitz([1/gamma(1-alpha)*B_l'*ws+d(0,1,alpha,sigma,tau)-2*(1-sigma)/h^2,(1-sigma)/h^2,zeros(M-3,1)'])*u(2:M,n-1)...
        +(1-sigma)/h^2*[u(1,n-1);zeros(M-3,1);u(M+1,n-1)]...
        +sigma/h^2*[u(1,n);zeros(M-3,1);u(M+1,n)];
    % 求解代数系统（由 k-2 层及第 k-1 层求第 k 层的值）
    u(2:M,n)=Coefficient_Matrix_A\Right_Term_B1;
    fprintf('进程:\t%d/%d\n',n,N+1)
end
% clc
fprintf('时间步长为 %f 空间步长为 %f 时快速 L2_1_sigma 插值逼近 %d\n',tau,h)
toc
%% 误差分析
U=zeros(size(u));
for m=1:M+1
    for n=1:N+1
        U(m,n)=u_exact(x(m),t(n),0,0);
    end
end
Error=max(max(abs(u-U)))

%% 绘图
figure
plot(t,u(7,:))
hold on
plot(t,U(7,:))
legend('数值解','精确解')

创作不易，此处仅展示了主程序代码部分, 绘图及误差分析等完整代码请查看专栏Matlab偏微分方程系列编程，感谢支持.

参考文献

孙志忠,高广花.分数阶微分方程的有限差分方法(第二版).北京：科学出版社,2021.

---

本人水平有限, 若有不妥之处, 恳请批评指正.

作者:图灵的猫

作者邮箱: turingscat@126.com