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

图灵猫

已于 2023-03-13 18:50:40 修改

阅读量3.2k

点赞数 9

分类专栏：偏微分方程

于 2022-06-20 10:27:02 首次发布

本文链接：https://blog.csdn.net/qq_42818403/article/details/125367493

版权

matlab 快速算法分数阶

偏微分方程专栏收录该内容

15 篇文章

订阅专栏

Caputo 分数阶一维问题基于 L1 逼近的快速差分方法

Caputo 分数阶一维问题基于 L1 逼近的空间二阶方法(附Matlab程序)

文章目录

Caputo 分数阶一维问题基于 L1 逼近的快速差分方法

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

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

本文给出求解上述时间分数阶慢扩散方程初边值问题基于 L1 逼近的快速差分方法.

差分格式的建立

在结点 $\left(x_{i}, t_{n}\right)$ 处考虑微分方程, 得到
$}_{0}^{C} D_{t}^{\alpha} u\left(x_{i}, t_{n}\right)=u_{x x}\left(x_{i}, t_{n}\right)+f_{i}^{n}, \quad 1 \leqslant i \leqslant M-1, \quad 1 \leqslant n \leqslant N .$
对上式中时间分数阶导数应用快速的 L1 公式离散、空间二阶导数应用二阶中心差商离散, 可得
$\left\{\begin{array}{c} \frac{1}{\Gamma(1-\alpha)}\left[\sum\limits_{l=1}^{N_{\exp }} \omega_{l} F_{l, i}^{n}+\hat{a}_{0}^{(\alpha)}\left(U_{i}^{n}-U_{i}^{n-1}\right)\right]=\delta_{x}^{2} U_{i}^{n}+f_{i}^{n}+\left(r_{4}\right)_{i}^{n}, \\ 1 \leqslant i \leqslant M-1, \quad 1 \leqslant n \leqslant N, \\ F_{l, i}^{1}=0, \quad F_{l, i}^{n}=\mathrm{e}^{-s_{l} \tau} F_{l, i}^{n-1}+B_{l}\left(U_{i}^{n-1}-U_{i}^{n-2}\right), \\ 1 \leqslant l \leqslant N_{\exp }, 1 \leqslant i \leqslant M-1,2 \leqslant n \leqslant N, \end{array}\right.$

且存在正常数 $c_{4}$ 使得
$\left|\left(r_{4}\right)_{i}^{n}\right| \leqslant c_{4}\left(\tau^{2-\alpha}+h^{2}+\epsilon\right), \quad 1 \leqslant i \leqslant M-1,1 \leqslant n \leqslant N .$
注意到初边值条件, 有
$\begin{cases}U_{i}^{0}=\varphi\left(x_{i}\right), & 1 \leqslant i \leqslant M-1, \\ U_{0}^{n}=\mu\left(t_{n}\right), & U_{M}^{n}=\nu\left(t_{n}\right), \quad 0 \leqslant n \leqslant N .\end{cases}$
在上式中略去小量项 $\left(r_{4}\right)_{i}^{n}$ , 并用数值解 $u_{i}^{n}$ 代替精确解 $U_{i}^{n}$ , 可得求解问题的如下差分格式:
$\left\{\begin{array}{l} \frac{1}{\Gamma(1-\alpha)}\left[\sum\limits_{l=1}^{N_{\exp }} \omega_{l} F_{l, i}^{n}+\hat{a}_{0}^{(\alpha)}\left(u_{i}^{n}-u_{i}^{n-1}\right)\right]=\delta_{x}^{2} u_{i}^{n}+f_{i}^{n} \\ 1 \leqslant i \leqslant M-1,1 \leqslant n \leqslant N \\ F_{l, i}^{1}=0, \quad F_{l, i}^{n}=\mathrm{e}^{-s_{l} \tau} F_{l, i}^{n-1}+B_{l}\left(u_{i}^{n-1}-u_{i}^{n-2}\right) \\ 1 \leqslant l \leqslant N_{\exp }, 1 \leqslant i \leqslant M-1,2 \leqslant n \leqslant N \\ u_{i}^{0}=\varphi\left(x_{i}\right), \quad 1 \leqslant i \leqslant M-1, \\ u_{0}^{n}=\mu\left(t_{n}\right), \quad u_{M}^{n}=\nu\left(t_{n}\right), \quad 0 \leqslant n \leqslant N \end{array}\right.$

当 $\left\{u_{i}^{k} \mid 0 \leqslant i \leqslant M, 0 \leqslant k \leqslant n-1\right\}$ 已知时, 得到 $\left\{F_{l, i}^{n} \mid 1 \leqslant l \leqslant N_{\exp }, 1 \leqslant i \leqslant M-1\right\}$ , 并将其代入差分格式的左端, 得到关于 $\left\{u_{i}^{n} \mid 0 \leqslant i \leqslant M\right\}$ 的线性方程组, 其运算量为 $O\left(M N_{\exp }\right)$ . 用追赶法解三对角方程组的运算量为 $O (M)$ . 因而此快速算法的总运算量为 $O\left(M N N_{\exp }\right)$ .

由传统的基于 L1 逼近得到关于 $\left\{u_{i}^{n} \mid 0 \leqslant i \leqslant M\right\}$ 的线性方程组, 其运算量为 $O (M n)$ . 其算法的总运算量为 $O\left(M \sum\limits_{n=1}^{N} n\right)=O\left(M N^{2}\right)$ .

当 $\gg N_{\exp }$ 时, $O\left(M N^{2}\right) \gg O\left(M N N_{\exp }\right)$ . 我们将上述算法称为基于 L1 逼近的快速算法或快速的 L1 差分格式.

分数阶微分方程数值算例

数值算例

考虑问题:
$}_{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}.$
则有
$\left\{\begin{array}{l}u_{x}=t^{2} \cos x \\ u_{x x}=-t^{2} \sin x\end{array} \quad\left\{\begin{array}{l}u_{t}=2 t \sin x \\ u_{t t}=2 \sin x .\end{array}\right.\right.$

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

$\begin{aligned} { }_{0}^{c} D_{t}^{\alpha} u(x, t)&=\frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} \frac{2s \sin x}{(t-s)^{\alpha}} d s\\ &=\frac{-2 \sin x}{\Gamma(1-\alpha)} \int_{0}^{t}(t-s-t)(t-s)^{-\alpha} d s\\ &=\frac{-2 \sin x}{\Gamma(1-\alpha)} \int_{0}^{t}\left[(t-s)^{1-\alpha}-t(t-s)^{-\alpha}\right] d s\\ &=\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[\left.\frac{1}{2-\alpha}(t-s)^{2-\alpha}\right|_{t} ^{0}+\left.\frac{t}{1-\alpha}(t-s)^{1-\alpha}\right|_{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{aligned}$

由于 $}_{0}^{C} D_{t}^{\alpha} u(x, t)=\frac{\partial^{2} u}{\partial x^{2}}(x, t)+f(x, t)$
于是 $\begin{aligned} 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{aligned}$

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

$\left\{\begin{array}{l} { }_{0}^{C} D_{t}^{\alpha} u(x, t)=u_{x x}(x, t)+\frac{2 \sin x}{\Gamma(3-\alpha)} t^{2-\alpha}+t^{2} \sin x, \quad x \in(0, L), t \in(0, T] \\ u(x, 0)=0, \quad x \in(0, L) \\ u(0, t)=0, \quad u(L, t)=t^2\sin L, \quad t \in[0, T] \end{array}\right.$

代数系统

由差分格式:
$\frac{1}{\Gamma(1-\alpha)}\left[\sum\limits_{l=1}^{N_{\exp }} \omega_{l} F_{l, i}^{n}+\hat{a}_{0}^{(\alpha)}\left(u_{i}^{n}-u_{i}^{n-1}\right)\right]=\delta_{x}^{2} u_{i}^{n}+f_{i}^{n}$
得到如下代数系统:
$\mathbf{A}\left(\begin{array}{c}u_{1}^{n} \\ u_{2}^{n} \\ \vdots \\ u_{M-2}^{n} \\ u_{M-1}^{n}\end{array}\right)=\mathbf{b}$
其中
$\mathbf{A}=\left(\begin{array}{cccc}\frac{1}{\Gamma(1-\alpha)} \hat{a}_{0}^{(\alpha)}+\frac{2}{h^{2}}&-\frac{1}{h^{2}} & & \\ -\frac{1}{h^{2}} & \frac{1}{\Gamma(1-\alpha)} \hat{a}_{0}^{(\alpha)}+\frac{2}{h^{2}}& -\frac{1}{h^{2}} & \\ & \ddots & \ddots & \\ & & \frac{1}{\Gamma(1-\alpha)} \hat{a}_{0}^{(\alpha)}+\frac{2}{h^{2}}& -\frac{1}{h^{2}}\end{array}\right)$

$\mathbf{b}=\frac{1}{h^{2}} \cdot\left(\begin{array}{c}u_{0}^{n} \\ 0 \\ \vdots \\ 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} \\ \vdots \\ u_{M-2}^{n-1}\\ u_{M-1}^{n-1}\end{array}\right)-\frac{1}{\Gamma(1-\alpha)} \sum\limits_{l=1}^{N_{exp}} \omega_{l} F_{l, i}^{n}+f_{i}^{n}$

数值结果

时间方向参数选取:
FL1参数设置-时间方向-

时间方向数值结果:
FL1数值结果-时间方向-

空间方向参数选取:
FL1参数设置-空间方向-

空间方向数值结果:
FL1数值结果-空间方向-

程序运行时间:

时间步长为 0.000100 空间步长为 0.012500 时快速 L1 插值逼近历时 4.104740 秒.

误差Error =6.9887e-07

而相同的网格剖分及精度下, 传统的基于 L1 逼近需历时 122.812862 秒.

误差Error =6.9845e-07

源代码

主程序:

clc,clear
tic
%% 初始化定解问题
alpha=0.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([2/h^2+1/gamma(1-alpha)*a_hat(0,alpha,tau),-1/h^2,zeros(M-3,1)']);
% 右端 Right_Term_B
Right_Term_B=1/gamma(1-alpha)*a_hat(0,alpha,tau)*u(2:M,n-1)...
    +f_source(x(2:M),t(n),alpha)'...
    +1/h^2.*[u(1,n);zeros(M-3,1);u(M+1,n)]...
    -1/gamma(1-alpha)*F_l*ws;
% 求解迭代系统（由0层求第1层的值）
u(2:M,n)=Coefficient_Matrix_A\Right_Term_B;
fprintf('进程:\t%d/%d\n',n,N+1)

%% 三层格式
B_l=1./(-xs*tau).*(exp(1).^(-2*xs*tau)-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)))*B_l';
    % 系数矩阵 Coefficient_Matrix_A
    Coefficient_Matrix_A=toeplitz([2/h^2+1/gamma(1-alpha)*a_hat(0,alpha,tau),-1/h^2,zeros(M-3,1)']);
    % 右端 Right_Term_B1
    Right_Term_B1=1/gamma(1-alpha)*a_hat(0,alpha,tau)*u(2:M,n-1)...
        +f_source(x(2:M),t(n),alpha)'...
        +1/h^2.*[u(1,n);zeros(M-3,1);u(M+1,n)]...
        -1/gamma(1-alpha)*F_l*ws;
    % 求解代数系统（由 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 时快速 L1 插值逼近 %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('数值解','精确解')