Caputo 分数阶导数快速的 H2N2 插值逼近(附Matlab代码)

Caputo 分数阶导数快速的 H2N2 插值逼近

Caputo 分数阶导数的 H2N2 插值逼近 (附Matlab程序)

H2N2 插值逼近格式推导

众所周知, 分数阶微分方程由于其存在全局依赖(非局部性)加大了求解困难, 于是有学者给出了一系列的快速算法,此处我们讨论 Caputo 分数阶导数基于 H2N2 逼近的快速算法, 其中分数阶导数的阶 $\gamma \in (1,2)$.

本文中, 为简单起见, 省去了 $N_{\exp }^{(\gamma -1)},s_l^{(\gamma -1)},{\omega }_l^{(\gamma -1)}$ 中的上标 $(\gamma -1)$. 应用引理 1.7.1 (见文末参考文献) 可得
$$\begin{array}{rlr}& {{}_0^C{D}_t^{\gamma }f(t)|}_{t=t_{n-\frac{1}{2}}}& \\ =& \frac{1}{\Gamma (2-\gamma )}\left[{\int }_{t_0}^{t_{\frac{1}{2}}}{f}^{\prime \prime }(t){\left(t_{n-\frac{1}{2}}-t\right)}^{1-\gamma }\mathrm{d}t+\sum _{k=1}^{n-1}{\int }_{t_{k-\frac{1}{2}}}^{t_{k+\frac{1}{2}}}{f}^{\prime \prime }(t){\left(t_{n-\frac{1}{2}}-t\right)}^{1-\gamma }\mathrm{d}t\right]& \\ \approx & \frac{1}{\Gamma (2-\gamma )}[{\int }_{t_0}^{t_{\frac{1}{2}}}{H}_{2,0}^{\prime \prime }(t)\sum _{l=1}^{N_{\text{exp }}}{\omega }_l{\mathrm{e}}^{-s_l\left(t_{n-\frac{1}{2}}-t\right)}\mathrm{d}t& +\sum _{k=1}^{n-2}{\int }_{t_{k-\frac{1}{2}}}^{t_{k+\frac{1}{2}}}{N}_{2,k}^{\prime \prime }(t)\sum _{l=1}^{N_{\text{exp }}}{\omega }_l{\mathrm{e}}^{-s_l\left(t_{n-\frac{1}{2}}-t\right)}\mathrm{d}t+{\int }_{t_{n-\frac{3}{2}}}^{t_{n-\frac{1}{2}}}{N}_{2,n-1}^{\prime \prime }(t){\left(t_{n-\frac{1}{2}}-t\right)}^{1-\gamma }\mathrm{d}t]\\ =& \frac{1}{\Gamma (2-\gamma )}\left[\sum _{l=1}^{N_{\text{exp }}}{\omega }_l{F}_l^n+{\delta }_t^2{f}^{n-1}{\int }_{t_{n-\frac{3}{2}}}^{t_{n-\frac{1}{2}}^2}{\left(t_{n-\frac{1}{2}}-t\right)}^{1-\gamma }\mathrm{d}t\right]& \\ \equiv & {}^{\mathcal{F}}{\hat{D}}^{\gamma }f\left(t_{n-\frac{1}{2}}\right),2⩽n⩽N,& \end{array}$$
其中
$$\begin{array}{rl}{F}_l^n={\int }_{t_0}^{t_{\frac{1}{2}}}{H}_{2,0}^{\prime \prime }(t){\mathrm{e}}^{-s_l\left(t_{n-\frac{1}{2}}-t\right)}\mathrm{d}t+& \sum _{k=1}^{n-2}{\int }_{t_{k-\frac{1}{2}}}^{t_{k+\frac{1}{2}}}{N}_{2,k}^{\prime \prime }(t){\mathrm{e}}^{-s_l\left(t_{n-\frac{1}{2}}-t\right)}\mathrm{d}t\\ & 1⩽l⩽N_{\exp },2⩽n⩽N\end{array}$$

利用如下递推算法计算 $F_l^n$ :
$F_l^n={\int }_{t_0}^{t_{\frac{1}{2}}}{H}_{2,0}^{\prime \prime }(t){\mathrm{e}}^{-s_l\left(t_{n-\frac{1}{2}}-t\right)}\mathrm{d}t+{\sum }_{k=1}^{n-2}{\int }_{t_{k-\frac{1}{2}}}^{t_{k+\frac{1}{2}}}{N}_{2,k}^{\prime \prime }(t){\mathrm{e}}^{-s_l\left(t_{n-\frac{1}{2}}-t\right)}\mathrm{d}t$
$={\mathrm{e}}^{-s_l\tau }\left[{\int }_{t_0}^{t_1}{H}_{2,0}^{\prime \prime }(t){\mathrm{e}}^{-s_l\left(t_{n-\frac{3}{2}}-t\right)}\mathrm{d}t+{\sum }_{k=1}^{n-3}{\int }_{t_{k-\frac{1}{2}}}^{t_{k+\frac{1}{2}}}{N}_{2,k}^{\prime \prime }(t){\mathrm{e}}^{-s_l\left(t_{n-\frac{3}{2}}-t\right)}\mathrm{d}t\right]$
$+{\int }_{t_{n-\frac{5}{2}}}^{t_{n-\frac{3}{2}}}{N}_{2,n-2}^{\prime \prime }(t){\mathrm{e}}^{-s_l\left(t_{n-\frac{1}{2}}-t\right)}\mathrm{d}t$
$={\mathrm{e}}^{-s_l\tau }{F}_l^{n-1}+{\delta }_t^2{f}^{n-2}{\int }_{t_{n-\frac{5}{2}}}^{t_{n-\frac{3}{2}}}{\mathrm{e}}^{-s_l\left(t_{n-\frac{1}{2}}-t\right)}\mathrm{d}t$
$={\mathrm{e}}^{-s_l\tau }{F}_l^{n-1}+B_l\left({\delta }_t{f}^{n-\frac{3}{2}}-{\delta }_t{f}^{n-\frac{5}{2}}\right),3⩽n⩽N$,
其中
$$F_l^2={\int }_{t_0}^{t_{\frac{1}{2}}}{H}_{2,0}^{\prime \prime }(t){\mathrm{e}}^{-s_l\left(t_{\frac{3}{2}}-t\right)}\mathrm{d}t=\frac{2}{s_l\tau }\left({\mathrm{e}}^{-s_l\tau }-{\mathrm{e}}^{-\frac{3}{2}{s}_l\tau }\right)\left({\delta }_t{f}^{\frac{1}{2}}-f^{\prime }\left(t_0\right)\right)$$
$$B_l=\frac{1}{\tau }{\int }_{t_{n-\frac{5}{2}}}^{t_{n-\frac{3}{2}}}{\mathrm{e}}^{-s_l\left(t_{n-\frac{1}{2}}-t\right)}\mathrm{d}t=\frac{1}{s_l\tau }\left({\mathrm{e}}^{-s_l\tau }-{\mathrm{e}}^{-2s_l\tau }\right)$$

综上得到计算 ${}_0^C{D}_t^{\gamma }f(t){\mid }_{t=t_{n-\frac{1}{2}}}$ 的如下算法:
$$\{\begin{array}{l}{}^F{\hat{D}}^{\gamma }f\left(t_{n-\frac{1}{2}}\right)=\frac{1}{\Gamma (2-\gamma )}\left[\sum _{l=1}^{N_{\exp }}{\omega }_l{F}_l^n+\frac{{\tau }^{2-\gamma }}{2-\gamma }{\delta }_t^2{f}^{n-1}\right],2⩽n⩽N,\\ F_l^2=\frac{2}{\tau }{\int }_{t_0}^{t_{\frac{1}{2}}}{\mathrm{e}}^{-s_l\left(t_{\frac{3}{2}}-t\right)}\mathrm{d}t\left({\delta }_t{f}^{\frac{1}{2}}-f^{\prime }\left(t_0\right)\right),1⩽l⩽N_{\exp }\\ F_l^n={\mathrm{e}}^{-s_l\tau }{F}_l^{n-1}+B_l\left({\delta }_t{f}^{n-\frac{3}{2}}-{\delta }_t{f}^{n-\frac{5}{2}}\right),1⩽l⩽N_{\exp },3⩽n⩽N.\end{array}$$
我们称上述公式为快速的 H2N2 逼近或快速的 H2N2 公式. 应用快速的 H2N2 公式可以对时间分数阶波方程建立 $3-\gamma$ 阶时间精度的差分格式.

分数阶微分方程数值算例

数值算例

求解初值问题
$$\{\begin{array}{l}{}_0^C{\mathbf{D}}_t^{\gamma }u(t)=f(t),0<t⩽T\\ u(0)=0,u^{\prime }(0)=0\end{array}$$
其中 $\gamma \in (1,2)$, $f(t)=\frac{6}{\Gamma (2-\gamma )}(\frac{1}{2-\gamma }-\frac{1}{3-\gamma })t^{3-\gamma }.$

该方程精确解为 $u(t)=t^3.$

迭代格式

将数值算例结合 H2N2 逼近公式可得到如下迭代格式.

启动层(两层格式)

由于启动层不存在明显全局依赖, 故直接采用如下的传统 H2N2 格式求解即可.
$$\frac{1}{\Gamma (2-\gamma )}\left[{\hat{b}}_0^{(1,\gamma )}\cdot \frac{u^1-u^0}{\tau }-{\hat{b}}_0^{(1,\gamma )}{u}^{\prime }\left(t_0\right)\right]=f\left(t_1\right)$$
$$\widehat{b_0^{(1,\gamma )}}\cdot \frac{u^1-u^0}{\tau }-{\hat{b}}_0^{(1,\gamma )}{u}^{\prime }\left(t_0\right)=f\left(t_1\right)\cdot \Gamma (2-r)$$
得到启动层:
$$u\left(t_1\right)=\frac{\tau \cdot \Gamma (2-r)}{{\hat{b}}_0^{(1,\gamma )}}\cdot f\left(t_1\right)+\tau {\hat{b}}_0^{(1,\gamma )}{u}^{\prime }\left(t_0\right)+u\left(t_0\right).$$

三层格式

$$\frac{1}{\Gamma (2-\gamma )}[\sum _{l=1}^{N_{\text{exp }}}{\omega }_l{F}_l^n+\frac{{\tau }^{2-\gamma }}{2-\gamma }\left(\frac{1}{{\tau }^2}{u}^n-\frac{2}{{\tau }^2}{u}^{n-1}+\frac{1}{{\tau }^2}{u}^{n-2}\right)=f\left(t_{n-\frac{1}{2}}\right)$$

得到三层迭代格式:

$$u\left(t_n\right)=\Gamma (3-\gamma )\cdot {\tau }^{\gamma }\cdot f\left(t_{n-\frac{1}{2}}\right)-(2-\gamma ){\tau }^{\gamma }\sum _{i=1}^{N_{exp}}{\omega }_l{F}_l^n+2u\left(t_{n-1}\right)-u\left(t_{n-2}\right)$$

数值结果

参数选取:

数值结果:

程序运行时间:

时间步长为 0.000100 时快速 H2N2 插值逼近历时 0.751304 秒.

而相同的网格剖分及精度下, 传统的 H2N2 插值逼近需历时 57.789599 秒.

源代码

主程序:

clc,clear
tic
%% 初始化定解问题
alpha=1.2;
epsilon=1e-8;
tau_hat=1e-6;
tau=1/10000;           %   @tau 时间步长            
T_a=0;               %   @T 时间求解区域
T_b=1;
N=(T_b-T_a)/tau;     %   @N t被分割的区间数
t=T_a:tau:T_b;       %   @t 时间向量
u=zeros(1,N+1);      %   @u 初始化数值解
u(1)=f_ic(t(1),0);   %   定义初值条件
[xs,ws,nexp] = sumofexpappr2new(alpha-1,epsilon,tau_hat,T_b);
%% 两层格式（启动层）
n=2;
m=n-1;
% 求解迭代系统（由0层求第1层的值）
u(n)=tau*gamma(2-alpha)/b_hat(0,m,alpha,tau)*f_source((t(2)+t(1))/2,alpha)...
    +tau*b_hat(0,m,alpha,tau)*f_ic(t(1),1)...
    +u(n-1);
 
fprintf('进程:\t%d/%d\n',n,N+1)

%% 三层格式
F_l=2./(xs*tau).*(exp(1).^(-xs*tau)-exp(1).^(-3/2*xs*tau)).*(1/tau*(u(2)-u(1))-f_ic(t(1),1));
B_l=1./(xs*tau).*(exp(1).^(-xs*tau)-exp(1).^(-2*xs*tau));
for n=3:N+1
    % 求解迭代系统（由 k-2 层及第 k-1 层求第 k 层的值）
    u(n)=gamma(3-alpha)*tau^alpha*f_source((t(n)+t(n-1))/2,alpha)...
        -(2-alpha)*tau^alpha*(ws'*F_l)...
        +2*u(n-1)-u(n-2);
    F_l=exp(1).^(-xs*tau).*F_l+B_l.*(1/tau*(u(n)-2*u(n-1)+u(n-2)));
    fprintf('进程:\t%d/%d\n',n,N+1)
end
clc
fprintf('时间步长为 %f 时快速 H2N2 插值逼近 %d\n',tau)
toc
%% 误差分析
U=zeros(size(u));
for n=1:N+1
        U(n)=u_exact(t(n),0);
end
Error=max(max(abs(u-U)))

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

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

参考文献

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

---

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

作者:图灵的猫

作者邮箱: turingscat@126.com