Caputo 分数阶导数的 H2N2 插值逼近
Caputo 分数阶导数快速的 H2N2 插值逼近(附Matlab代码)
H2N2 插值逼近格式推导
本文介绍 H2N2 插值逼近方法, 考虑
0
C
D
t
γ
f
(
t
)
∣
t
=
t
n
−
1
2
=
1
Γ
(
2
−
γ
)
∫
t
0
t
n
−
1
2
f
′
′
(
t
)
(
t
n
−
1
2
−
t
)
1
−
γ
d
t
,
1
⩽
n
⩽
N
\left.{ }_{0}^{C} D_{t}^{\gamma} f(t)\right|_{t=t_{n-\frac{1}{2}}}=\frac{1}{\Gamma(2-\gamma)} \int_{t_{0}}^{t_{n-\frac{1}{2}}} f^{\prime \prime}(t)\left(t_{n-\frac{1}{2}}-t\right)^{1-\gamma} \mathrm{d} t, \quad 1 \leqslant n \leqslant N
0CDtγf(t)
t=tn−21=Γ(2−γ)1∫t0tn−21f′′(t)(tn−21−t)1−γdt,1⩽n⩽N
的近似计算. 将其写成多个小区间上和的形式可得
0
C
D
t
γ
f
(
t
)
∣
t
=
t
n
−
1
2
=
1
Γ
(
2
−
γ
)
[
∫
t
0
t
1
f
′
′
(
t
)
(
t
n
−
1
2
−
t
)
1
−
γ
d
t
+
∑
k
=
1
n
−
1
∫
t
k
−
1
2
t
k
+
1
2
f
′
′
(
t
)
(
t
n
−
1
2
−
t
)
1
−
γ
d
t
]
.
\begin{aligned} &\left.{ }_{0}^{C} D_{t}^{\gamma} f(t)\right|_{t=t_{n-\frac{1}{2}}} \\ =& \frac{1}{\Gamma(2-\gamma)}\left[\int_{t_{0}}^{t_{1}} 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] . \end{aligned}
=0CDtγf(t)
t=tn−21Γ(2−γ)1
∫t0t1f′′(t)(tn−21−t)1−γdt+k=1∑n−1∫tk−21tk+21f′′(t)(tn−21−t)1−γdt
.
利用
(
t
0
,
f
(
t
0
)
)
,
(
t
0
,
f
′
(
t
0
)
)
,
(
t
1
,
f
(
t
1
)
)
\left(t_{0}, f\left(t_{0}\right)\right),\left(t_{0}, f^{\prime}\left(t_{0}\right)\right),\left(t_{1}, f\left(t_{1}\right)\right)
(t0,f(t0)),(t0,f′(t0)),(t1,f(t1)) 作
f
(
t
)
f(t)
f(t) 的二次 Hermite 插值多项式得到
H
2
,
0
(
t
)
=
f
(
t
0
)
+
f
′
(
t
0
)
(
t
−
t
0
)
+
1
τ
(
δ
t
f
1
2
−
f
′
(
t
0
)
)
(
t
−
t
0
)
2
.
H_{2,0}(t)=f\left(t_{0}\right)+f^{\prime}\left(t_{0}\right)\left(t-t_{0}\right)+\frac{1}{\tau}\left(\delta_{t} f^{\frac{1}{2}}-f^{\prime}\left(t_{0}\right)\right)\left(t-t_{0}\right)^{2} .
H2,0(t)=f(t0)+f′(t0)(t−t0)+τ1(δtf21−f′(t0))(t−t0)2.
利用三点
(
t
k
−
1
,
f
(
t
k
−
1
)
)
,
(
t
k
,
f
(
t
k
)
)
,
(
t
k
+
1
,
f
(
t
k
+
1
)
)
\left(t_{k-1}, f\left(t_{k-1}\right)\right),\left(t_{k}, f\left(t_{k}\right)\right),\left(t_{k+1}, f\left(t_{k+1}\right)\right)
(tk−1,f(tk−1)),(tk,f(tk)),(tk+1,f(tk+1)) 作
f
(
t
)
f(t)
f(t) 的二次 Newton 插值多项式
N
2
,
k
(
t
)
=
f
(
t
k
−
1
)
+
(
δ
t
f
k
−
1
2
)
(
t
−
t
k
−
1
)
+
1
2
(
δ
t
2
f
k
)
(
t
−
t
k
−
1
)
(
t
−
t
k
)
,
N_{2, k}(t)=f\left(t_{k-1}\right)+\left(\delta_{t} f^{k-\frac{1}{2}}\right)\left(t-t_{k-1}\right)+\frac{1}{2}\left(\delta_{t}^{2} f^{k}\right)\left(t-t_{k-1}\right)\left(t-t_{k}\right),
N2,k(t)=f(tk−1)+(δtfk−21)(t−tk−1)+21(δt2fk)(t−tk−1)(t−tk),
记
b
^
k
(
n
,
γ
)
=
{
τ
1
−
γ
2
−
γ
[
(
k
+
1
)
2
−
γ
−
k
2
−
γ
]
,
0
⩽
k
⩽
n
−
2
2
τ
1
−
γ
2
−
γ
[
(
n
−
1
2
)
2
−
γ
−
(
n
−
1
)
2
−
γ
]
,
k
=
n
−
1.
\hat{b}_{k}^{(n, \gamma)}= \begin{cases}\frac{\tau^{1-\gamma}}{2-\gamma}\left[(k+1)^{2-\gamma}-k^{2-\gamma}\right], & 0 \leqslant k \leqslant n-2 \\ \frac{2 \tau^{1-\gamma}}{2-\gamma}\left[\left(n-\frac{1}{2}\right)^{2-\gamma}-(n-1)^{2-\gamma}\right], & k=n-1 .\end{cases}
b^k(n,γ)={2−γτ1−γ[(k+1)2−γ−k2−γ],2−γ2τ1−γ[(n−21)2−γ−(n−1)2−γ],0⩽k⩽n−2k=n−1.
计算可得
2
τ
∫
t
0
t
1
(
t
n
−
1
1
−
t
)
1
−
γ
d
t
=
2
τ
⋅
1
2
−
γ
[
(
t
n
−
1
2
−
t
0
)
2
−
γ
−
(
t
n
−
1
2
−
t
1
2
)
2
−
γ
]
=
2
τ
1
−
γ
2
−
γ
[
(
n
−
1
2
)
2
−
γ
−
(
n
−
1
)
2
−
γ
]
=
b
^
n
−
1
(
n
,
γ
)
\begin{aligned} & \frac{2}{\tau} \int_{t_{0}}^{t_{1}}\left(t_{n-\frac{1}{1}}-t\right)^{1-\gamma} \mathrm{d} t \\ =& \frac{2}{\tau} \cdot \frac{1}{2-\gamma}\left[\left(t_{n-\frac{1}{2}}-t_{0}\right)^{2-\gamma}-\left(t_{n-\frac{1}{2}}-t_{\frac{1}{2}}\right)^{2-\gamma}\right] \\ =& \frac{2 \tau^{1-\gamma}}{2-\gamma}\left[\left(n-\frac{1}{2}\right)^{2-\gamma}-(n-1)^{2-\gamma}\right] \\ =& \hat{b}_{n-1}^{(n, \gamma)} \end{aligned}
===τ2∫t0t1(tn−11−t)1−γdtτ2⋅2−γ1[(tn−21−t0)2−γ−(tn−21−t21)2−γ]2−γ2τ1−γ[(n−21)2−γ−(n−1)2−γ]b^n−1(n,γ)
和
1
τ
∫
t
k
−
1
2
t
k
+
1
2
(
t
n
−
1
2
−
t
)
1
−
γ
d
t
=
1
τ
⋅
1
2
−
γ
[
(
t
n
−
1
2
−
t
k
−
1
2
)
2
−
γ
−
(
t
n
−
1
2
−
t
k
+
1
2
)
2
−
γ
]
=
τ
1
−
γ
2
−
γ
[
(
n
−
k
)
2
−
γ
−
(
n
−
k
−
1
)
2
−
γ
]
=
b
^
n
−
k
−
1
(
n
,
γ
)
,
1
⩽
k
⩽
n
−
1.
\begin{aligned} & \frac{1}{\tau} \int_{t_{k-\frac{1}{2}}}^{t_{k+\frac{1}{2}}}\left(t_{n-\frac{1}{2}}-t\right)^{1-\gamma} \mathrm{d} t \\ =& \frac{1}{\tau} \cdot \frac{1}{2-\gamma}\left[\left(t_{n-\frac{1}{2}}-t_{k-\frac{1}{2}}\right)^{2-\gamma}-\left(t_{n-\frac{1}{2}}-t_{k+\frac{1}{2}}\right)^{2-\gamma}\right] \\ =& \frac{\tau^{1-\gamma}}{2-\gamma}\left[(n-k)^{2-\gamma}-(n-k-1)^{2-\gamma}\right] \\ =& \hat{b}_{n-k-1}^{(n, \gamma)}, \quad 1 \leqslant k \leqslant n-1 . \end{aligned}
===τ1∫tk−21tk+21(tn−21−t)1−γdtτ1⋅2−γ1[(tn−21−tk−21)2−γ−(tn−21−tk+21)2−γ]2−γτ1−γ[(n−k)2−γ−(n−k−1)2−γ]b^n−k−1(n,γ),1⩽k⩽n−1.
于是得到如下逼近公式
D
^
t
γ
f
(
t
n
−
1
2
)
=
1
Γ
(
2
−
γ
)
[
b
^
n
−
1
(
n
,
γ
)
⋅
(
δ
t
f
1
2
−
f
′
(
t
0
)
)
+
∑
k
=
1
n
−
1
b
^
n
−
k
−
1
(
n
,
γ
)
⋅
(
δ
t
f
k
+
1
2
−
δ
t
f
k
−
1
2
)
]
=
1
Γ
(
2
−
γ
)
[
b
^
0
(
n
,
γ
)
δ
t
f
n
−
1
2
−
∑
k
=
1
n
−
1
(
b
^
n
−
k
−
1
(
n
,
γ
)
−
b
^
n
−
k
(
n
,
γ
)
)
δ
t
f
k
−
1
2
−
b
^
n
−
1
(
n
,
γ
)
f
′
(
t
0
)
]
\begin{aligned} & \hat{D}_{t}^{\gamma} f\left(t_{n-\frac{1}{2}}\right) \\ =& \frac{1}{\Gamma(2-\gamma)}\left[\hat{b}_{n-1}^{(n, \gamma)} \cdot\left(\delta_{t} f^{\frac{1}{2}}-f^{\prime}\left(t_{0}\right)\right)+\sum_{k=1}^{n-1} \hat{b}_{n-k-1}^{(n, \gamma)} \cdot\left(\delta_{t} f^{k+\frac{1}{2}}-\delta_{t} f^{k-\frac{1}{2}}\right)\right] \\ =& \frac{1}{\Gamma(2-\gamma)}\left[\hat{b}_{0}^{(n, \gamma)} \delta_{t} f^{n-\frac{1}{2}}-\sum_{k=1}^{n-1}\left(\hat{b}_{n-k-1}^{(n, \gamma)}-\hat{b}_{n-k}^{(n, \gamma)}\right) \delta_{t} f^{k-\frac{1}{2}}-\hat{b}_{n-1}^{(n, \gamma)} f^{\prime}\left(t_{0}\right)\right] \end{aligned}
==D^tγf(tn−21)Γ(2−γ)1[b^n−1(n,γ)⋅(δtf21−f′(t0))+k=1∑n−1b^n−k−1(n,γ)⋅(δtfk+21−δtfk−21)]Γ(2−γ)1[b^0(n,γ)δtfn−21−k=1∑n−1(b^n−k−1(n,γ)−b^n−k(n,γ))δtfk−21−b^n−1(n,γ)f′(t0)]
我们称上述公式为 H2N2 逼近或 H2N2 公式. 应用 H2N2 公式可以对时间分数阶波方程建立
3
−
γ
3-\gamma
3−γ 阶时间精度的差分格式.
分数阶微分方程数值算例
数值算例
求解初值问题
{
0
C
D
t
γ
u
(
t
)
=
f
(
t
)
,
0
<
t
⩽
T
u
(
0
)
=
0
,
u
′
(
0
)
=
0
\left\{\begin{array}{l} { }_{0}^{C} \mathbf{D}_{t}^{\gamma} u(t)=f(t), \quad 0<t \leqslant T \\ u(0)=0, \quad u^{\prime}(0)=0 \end{array}\right.
{0CDtγu(t)=f(t),0<t⩽Tu(0)=0,u′(0)=0
其中
γ
∈
(
1
,
2
)
\gamma \in(1,2)
γ∈(1,2),
f
(
t
)
=
6
Γ
(
2
−
γ
)
(
1
2
−
γ
−
1
3
−
γ
)
t
3
−
γ
.
f(t)=\frac{6}{\Gamma(2-\gamma)}(\frac{1}{2-\gamma}-\frac{1}{3-\gamma})t^{3-\gamma}.
f(t)=Γ(2−γ)6(2−γ1−3−γ1)t3−γ.
该方程精确解为 u ( t ) = t 3 . u(t)=t^3. u(t)=t3.
迭代格式
将数值算例结合 H2N2 逼近公式可得到如下迭代格式.
启动层(两层格式)
1
Γ
(
2
−
γ
)
[
b
^
0
(
1
,
γ
)
⋅
u
1
−
u
0
τ
−
b
^
0
(
1
,
γ
)
u
′
(
t
0
)
]
=
f
(
t
1
)
\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)
Γ(2−γ)1[b^0(1,γ)⋅τu1−u0−b^0(1,γ)u′(t0)]=f(t1)
b
0
(
1
,
γ
)
^
⋅
u
1
−
u
0
τ
−
b
^
0
(
1
,
γ
)
u
′
(
t
0
)
=
f
(
t
1
)
⋅
Γ
(
2
−
r
)
\hat{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)
b0(1,γ)^⋅τu1−u0−b^0(1,γ)u′(t0)=f(t1)⋅Γ(2−r)
得到启动层:
u
(
t
1
)
=
τ
⋅
Γ
(
2
−
r
)
b
^
0
(
1
,
γ
)
⋅
f
(
t
1
)
+
τ
b
^
0
(
1
,
γ
)
u
′
(
t
0
)
+
u
(
t
0
)
.
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) .
u(t1)=b^0(1,γ)τ⋅Γ(2−r)⋅f(t1)+τb^0(1,γ)u′(t0)+u(t0).
三层格式
1 Γ ( 2 − γ ) [ b ^ 0 ( n , γ ) δ t u n − 1 2 − ∑ k = 1 n − 1 ( b ^ n − k − 1 ( n , γ ) − b ^ n − k ( n , γ ) ) δ t u k − 1 2 − b ^ n − 1 ( n , γ ) u ′ ( t 0 ) ] = f ( t n ) \frac{1}{\Gamma(2-\gamma)}\left[\hat{b}_{0}^{(n, \gamma)} \delta_{t} u^{n-\frac{1}{2}}-\sum_{k=1}^{n-1}\left(\hat{b}_{n-k-1}^{(n, \gamma)}-\hat{b}_{n-k}^{(n, \gamma)}\right) \delta_{t} u^{k-\frac{1}{2}}-\hat{b}_{n-1}^{(n, \gamma)} u^{\prime}\left(t_{0}\right)\right]=f\left(t_{n}\right) Γ(2−γ)1[b^0(n,γ)δtun−21−k=1∑n−1(b^n−k−1(n,γ)−b^n−k(n,γ))δtuk−21−b^n−1(n,γ)u′(t0)]=f(tn)
b ^ 0 ( n , γ ) ⋅ u n − u n − 1 τ − ∑ k = 1 n − 1 ( b ^ n , k − 1 ( n , γ ) − b ^ n − k ( n , γ ) ) u k − u k − 1 τ − b ^ n − 1 ( n , γ ) u ′ ( t 0 ) = Γ ( 2 − γ ) ⋅ f ( t n ) \hat{b}_{0}^{(n, \gamma)} \cdot \frac{u^{n}-u^{n-1}}{\tau}-\sum_{k=1}^{n-1}\left(\hat{b}_{n, k-1}^{(n, \gamma)}-\hat{b}_{n-k}^{(n, \gamma)}\right) \frac{u^{k}-u^{k-1}}{\tau}-\hat{b}_{n-1}^{(n, \gamma)} u^{\prime}\left(t_{0}\right)=\Gamma(2-\gamma) \cdot f\left(t_{n}\right) b^0(n,γ)⋅τun−un−1−k=1∑n−1(b^n,k−1(n,γ)−b^n−k(n,γ))τuk−uk−1−b^n−1(n,γ)u′(t0)=Γ(2−γ)⋅f(tn)
得到三层迭代格式:
u ( t n ) = τ Γ ( 2 − γ ) b ^ 0 ( n , γ ) ⋅ f ( t n ) + u ( t n − 1 ) + ∑ k = 1 n − 1 b ^ n − k − 1 ( n , γ ) − b ^ n − k ( n , γ ) b ^ 0 ( n , r ) ⋅ ( u ( t k ) − u ( t k − 1 ) ) + τ ⋅ b ^ n − 1 ( n , γ ) b ^ 0 ( n , γ ) f ′ ( t 0 ) u\left(t_{n}\right)=\frac{\tau\Gamma(2-\gamma)}{\hat{b}_{0}^{(n, \gamma)}} \cdot f\left(t_{n}\right)+u\left(t_{n-1}\right)+\sum_{k=1}^{n-1} \frac{\hat{b}_{n-k-1}^{(n, \gamma)} -\hat{b}_{n-k}^{(n, \gamma)}}{\hat{b}_{0}^{(n, r)}} \cdot\left(u\left(t_{k}\right)-u\left(t_{k-1}\right)\right)+\frac{\tau \cdot \hat{b}_{n-1}^{(n, \gamma)}}{\hat{b}_{0}^{(n, \gamma)}} f^{\prime}\left(t_{0}\right) u(tn)=b^0(n,γ)τΓ(2−γ)⋅f(tn)+u(tn−1)+k=1∑n−1b^0(n,r)b^n−k−1(n,γ)−b^n−k(n,γ)⋅(u(tk)−u(tk−1))+b^0(n,γ)τ⋅b^n−1(n,γ)f′(t0)
数值结果
参数选取:
数值结果:
程序运行时间:
时间步长为 0.000100 时 H2N2 插值逼近历时 57.789599
秒.
而相同的网格剖分及精度下, 快速的 H2N2 插值逼近仅需历时 0.751304
秒.
源代码
主程序:
clc,clear
tic
%% 初始化定解问题
alpha=1.2;
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); % 定义初值条件
%% 两层格式(启动层)
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)
%% 三层格式
for n=3:N+1
m=n-1;
% 生成求和项
S=0;
for k=1:m-1
S=S+(b_hat(m-k-1,m,alpha,tau)-b_hat(m-k,m,alpha,tau))...
/b_hat(0,m,alpha,tau)*(u(k+1)-u(k));
end
% 求解迭代系统(由 k-2 层及第 k-1 层求第 k 层的值)
u(n)=tau*gamma(2-alpha)/b_hat(0,m,alpha,tau)*f_source((t(n)+t(n-1))/2,alpha)...
+u(n-1)...
+S...
+tau*b_hat(m-1,m,alpha,tau)/b_hat(0,m,alpha,tau)*f_ic(t(1),1);
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