二阶导数差分表达式推导_二阶微分方程数学建模[通俗易懂]

二阶导数差分表达式推导_二阶微分方程数学建模[通俗易懂]本文介绍H2N2插值逼近方法,结合数值算例验证了收敛情况,并给出了Matlab程序源代码._caputo型分数阶导数

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=tn21
=
Γ(2γ)1t0tn21f′′(t)(tn21t)1γdt,1nN

的近似计算. 将其写成多个小区间上和的形式可得
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=tn21
Γ(2γ)1
t0t1f′′(t)(tn21t)1γdt+k=1n1tk21tk+21f′′(t)(tn21t)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)(tt0)+τ1(δtf21f(t0))(tt0)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) (tk1,f(tk1)),(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(tk1)+(δtfk21)(ttk1)+21(δt2fk)(ttk1)(ttk),


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γ[(n21)2γ(n1)2γ],0kn2k=n1.

计算可得
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} ===τ2t0t1(tn11t)1γdtτ22γ1[(tn21t0)2γ(tn21t21)2γ]2γ2τ1γ[(n21)2γ(n1)2γ]b^n1(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} ===τ1tk21tk+21(tn21t)1γdtτ12γ1[(tn21tk21)2γ(tn21tk+21)2γ]2γτ1γ[(nk)2γ(nk1)2γ]b^nk1(n,γ),1kn1.
于是得到如下逼近公式
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(tn21)Γ(2γ)1[b^n1(n,γ)(δtf21f(t0))+k=1n1b^nk1(n,γ)(δtfk+21δtfk21)]Γ(2γ)1[b^0(n,γ)δtfn21k=1n1(b^nk1(n,γ)b^nk(n,γ))δtfk21b^n1(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<tTu(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γ13γ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,γ)τu1u0b^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,γ)^τu1u0b^0(1,γ)u(t0)=f(t1)Γ(2r)
得到启动层:
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,γ)τΓ(2r)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,γ)δtun21k=1n1(b^nk1(n,γ)b^nk(n,γ))δtuk21b^n1(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,γ)τunun1k=1n1(b^n,k1(n,γ)b^nk(n,γ))τukuk1b^n1(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(tn1)+k=1n1b^0(n,r)b^nk1(n,γ)b^nk(n,γ)(u(tk)u(tk1))+b^0(n,γ)τb^n1(n,γ)f(t0)

数值结果

参数选取:
H2N2数值结果

数值结果:
H2N2数值结果图

程序运行时间:

时间步长为 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

今天的文章二阶导数差分表达式推导_二阶微分方程数学建模[通俗易懂]分享到此就结束了,感谢您的阅读。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
如需转载请保留出处:https://bianchenghao.cn/89358.html

(0)
编程小号编程小号

相关推荐

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注