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

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

Caputo 分数阶一维问题基于 L1 逼近的空间二阶方法(附Matlab程序)
文章目录
- Caputo 分数阶一维问题基于 L1 逼近的快速差分方法
- 差分格式的建立
  分数阶微分方程数值算例
  数值算例
  源项和初边值条件
  代数系统
  数值结果
  源代码
  
  参考文献
考虑如下时间分数阶慢扩散方程初边值问题
$\left\{$
$\begin{array}{l} _{0}^{C} D_{t}^{α} u (x, t) = u_{x x} (x, t) + f (x, t), x \in (0, L), t \in (0, T] \\ u (x, 0) = φ (x), x \in (0, L) \\ u (0, t) = μ (t), u (L, t) = ν (t), t \in [0, T] \end{array}$
\right. ⎩⎨⎧0CDtαu(x,t)=uxx(x,t)+f(x,t),x∈(0,L),t∈(0,T]u(x,0)=φ(x),x∈(0,L)u(0,t)=μ(t),u(L,t)=ν(t),t∈[0,T]
其中 $\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{matrix} \frac{1}{Γ (1 - α)} [\sum_{l = 1}^{N_{\exp}} ω_{l} F_{l, i}^{n} + {\hat{a}}_{0}^{(α)} (U_{i}^{n} - U_{i}^{n - 1})] = δ_{x}^{2} U_{i}^{n} + f_{i}^{n} + {(r_{4})}_{i}^{n}, \\ 1 ⩽ i ⩽ M - 1, 1 ⩽ n ⩽ N, \\ F_{l, i}^{1} = 0, F_{l, i}^{n} = e^{- s_{l} τ} F_{l, i}^{n - 1} + B_{l} (U_{i}^{n - 1} - U_{i}^{n - 2}), \\ 1 ⩽ l ⩽ N_{\exp}, 1 ⩽ i ⩽ M - 1, 2 ⩽ n ⩽ N, \end{matrix}$
\right. ⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧Γ(1−α)1[l=1∑NexpωlFl,in+a^0(α)(Uin−Uin−1)]=δx2Uin+fin+(r4)in,1⩽i⩽M−1,1⩽n⩽N,Fl,i1=0,Fl,in=e−slτFl,in−1+Bl(Uin−1−Uin−2),1⩽l⩽Nexp,1⩽i⩽M−1,2⩽n⩽N,

且存在正常数 $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} = φ (x_{i}), & 1 ⩽ i ⩽ M - 1, \\ U_{0}^{n} = μ (t_{n}), & U_{M}^{n} = ν (t_{n}), 0 ⩽ n ⩽ N . \end{cases}$
{Ui0=φ(xi),U0n=μ(tn),1⩽i⩽M−1,UMn=ν(tn),0⩽n⩽N.
在上式中略去小量项 $\left(r_{4}\right)_{i}^{n}$ , 并用数值解 $u_{i}^{n}$ 代替精确解 $U_{i}^{n}$ , 可得求解问题的如下差分格式:
$\begin{array}{l} \frac{1}{Γ (1 - α)} [\sum_{l = 1}^{N_{\exp}} ω_{l} F_{l, i}^{n} + {\hat{a}}_{0}^{(α)} (u_{i}^{n} - u_{i}^{n - 1})] = δ_{x}^{2} u_{i}^{n} + f_{i}^{n} \\ 1 ⩽ i ⩽ M - 1, 1 ⩽ n ⩽ N \\ F_{l, i}^{1} = 0, F_{l, i}^{n} = e^{- s_{l} τ} F_{l, i}^{n - 1} + B_{l} (u_{i}^{n - 1} - u_{i}^{n - 2}) \\ 1 ⩽ l ⩽ N_{\exp}, 1 ⩽ i ⩽ M - 1, 2 ⩽ n ⩽ N \\ u_{i}^{0} = φ (x_{i}), 1 ⩽ i ⩽ M - 1, \\ u_{0}^{n} = μ (t_{n}), u_{M}^{n} = ν (t_{n}), 0 ⩽ n ⩽ N \end{array}$

当 $\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. {ux=t2cosxuxx=−t2sinx{ut=2tsinxutt=2sinx.

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

$\begin{aligned} _{0}^{c} D_{t}^{α} u (x, t) & = \frac{1}{Γ (1 - α)} \int_{0}^{t} \frac{2 s \sin x}{(t - s)^{α}} d s \\ = \frac{- 2 \sin x}{Γ (1 - α)} \int_{0}^{t} (t - s - t) (t - s)^{- α} d s \\ = \frac{- 2 \sin x}{Γ (1 - α)} \int_{0}^{t} [(t - s)^{1 - α} - t (t - s)^{- α}] d s \\ = \frac{- 2 \sin x}{Γ (1 - α)} [\int_{0}^{t} - (t - s)^{1 - α} d (t - s) + \int_{0}^{1} t (t - s)^{- α} d (t - s)] \\ = \frac{- 2 \sin x}{Γ (1 - α)} [{\frac{1}{2 - α} (t - s)^{2 - α} |}_{t}^{0} + {\frac{t}{1 - α} (t - s)^{1 - α} |}_{0}^{t}] \\ = \frac{- 2 \sin x}{Γ (1 - α)} [\frac{t^{2 - α}}{2 - α} + \frac{t}{1 - α} \cdot (- t^{1 - α})] \\ = \frac{- 2 \sin x}{Γ (1 - α)} \cdot \frac{- t^{2 - α}}{(2 - α) (1 - α)} \\ = \frac{2 \sin x}{Γ (3 - α)} t^{2 - α} \end{aligned}$
0cDtαu(x,t)=Γ(1−α)1∫0t(t−s)α2ssinxds=Γ(1−α)−2sinx∫0t(t−s−t)(t−s)−αds=Γ(1−α)−2sinx∫0t[(t−s)1−α−t(t−s)−α]ds=Γ(1−α)−2sinx[∫0t−(t−s)1−αd(t−s)+∫01t(t−s)−αd(t−s)]=Γ(1−α)−2sinx[2−α1(t−s)2−α∣∣∣∣t0+1−αt(t−s)1−α∣∣∣∣0t]=Γ(1−α)−2sinx[2−αt2−α+1−αt⋅(−t1−α)]=Γ(1−α)−2sinx⋅(2−α)(1−α)−t2−α=Γ(3−α)2sinxt2−α

由于 $}_{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}^{α} u (x, t) - \frac{\partial^{2} u}{\partial x^{2}} (x, t) \\ = \frac{2 \sin x}{Γ (3 - α)} t^{2 - α} + t^{2} \sin x . \end{aligned}$
f(x,t)=0CDtαu(x,t)−∂x2∂2u(x,t)=Γ(3−α)2sinxt2−α+t2sinx.

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

$\left\{$
$\begin{array}{l} _{0}^{C} D_{t}^{α} u (x, t) = u_{x x} (x, t) + \frac{2 \sin x}{Γ (3 - α)} t^{2 - α} + 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}$
\right. ⎩⎨⎧0CDtαu(x,t)=uxx(x,t)+Γ(3−α)2sinxt2−α+t2sinx,x∈(0,L),t∈(0,T]u(x,0)=0,x∈(0,L)u(0,t)=0,u(L,t)=t2sinL,t∈[0,T]

代数系统

由差分格式:
$\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{matrix} u_{1}^{n} \\ u_{2}^{n} \\ ⋮ \\ u_{M - 2}^{n} \\ u_{M - 1}^{n} \end{matrix}$
\right)=\mathbf{b} A⎝⎜⎜⎜⎜⎜⎛u1nu2n⋮uM−2nuM−1n⎠⎟⎟⎟⎟⎟⎞=b
其中
$\begin{array}{cccc} \frac{1}{Γ (1 - α)} {\hat{a}}_{0}^{(α)} + \frac{2}{h^{2}} & - \frac{1}{h^{2}} \\ - \frac{1}{h^{2}} & \frac{1}{Γ (1 - α)} {\hat{a}}_{0}^{(α)} + \frac{2}{h^{2}} & - \frac{1}{h^{2}} \\ ⋱ & ⋱ \\ \frac{1}{Γ (1 - α)} {\hat{a}}_{0}^{(α)} + \frac{2}{h^{2}} & - \frac{1}{h^{2}} \end{array}$

$\mathbf{b}=\frac{1}{h^{2}} \cdot\left($
$\begin{matrix} u_{0}^{n} \\ 0 \\ ⋮ \\ 0 \\ u_{M}^{n} \end{matrix}$
\right)+\frac{\hat{a}_0^{(\alpha)}}{\Gamma(1-\alpha)} \cdot\left(
$\begin{matrix} u_{1}^{n - 1} \\ u_{2}^{n - 1} \\ ⋮ \\ u_{M - 2}^{n - 1} \\ u_{M - 1}^{n - 1} \end{matrix}$
\right)-\frac{1}{\Gamma(1-\alpha)} \sum\limits_{l=1}^{N_{exp}} \omega_{l} F_{l, i}^{n}+f_{i}^{n} b=h21⋅⎝⎜⎜⎜⎜⎜⎛u0n0⋮0uMn⎠⎟⎟⎟⎟⎟⎞+Γ(1−α)a^0(α)⋅⎝⎜⎜⎜⎜⎜⎛u1n−1u2n−1⋮uM−2n−1uM−1n−1⎠⎟⎟⎟⎟⎟⎞−Γ(1−α)1l=1∑NexpωlFl,in+fin

数值结果

时间方向参数选取:

时间方向数值结果:

空间方向参数选取:

空间方向数值结果:

程序运行时间:

时间步长为 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('数值解','精确解')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
```
此处由于字数限制, 仅展示了主程序代码部分, 若需要绘图及误差分析等完整代码, 请在评论区留下邮箱.

参考文献

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

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

作者:图灵的猫

作者邮箱: turingscat@126.com
相关阅读:
vector的模拟实现
 【每日刷题】Day22
第8章注意力机制与外部记忆
 零基础看深度学习python代码的基本方法
 [深度学习]基于C++和onnxruntime部署yolov10的onnx模型
 C语言第六课：函数（1）——分类、参数与调用
 Flink sql 实现 -connection-clickhouse的 source和 sink
import导入顺序杂乱的问题
 动态规划---最长序列问题详解
 Kotlin 的锁和多线程同步
原文地址：https://blog.csdn.net/qq_42818403/article/details/125367493

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

文章目录

差分格式的建立

分数阶微分方程数值算例

数值算例

源项和初边值条件

代数系统

数值结果

源代码

参考文献