四階時(shí)間分?jǐn)?shù)波方程的快速緊致差分方法

王 婉，張海湘，楊雪花

（湖南工業(yè)大學(xué) 理學(xué)院，湖南 株洲 412007）

1 研究背景

分?jǐn)?shù)微分方程現(xiàn)已被應(yīng)用于越來(lái)越多的領(lǐng)域，如黏彈性材料[1-2]、控制理論[3]、金融市場(chǎng)中的期權(quán)定價(jià)模型[4-5]等。到目前為止，科研工作者們提出了很多解決具有弱奇異性分?jǐn)?shù)階微分方程的方法，比如有限差分法[6-7]、有限元方法[8-9]、光譜方法[10]，傅里葉變換法[11-12]等。對(duì)于分?jǐn)?shù)階微分方程中的Caputo導(dǎo)數(shù)項(xiàng)，學(xué)者們采用不同方法近似處理。例如，Xu D.等[13]用L1離散Caputo導(dǎo)數(shù)項(xiàng)，用二階卷積求積公式近似Riemann-Liouville積分項(xiàng)，構(gòu)造了緊差分格式，不僅給出了收斂性和穩(wěn)定性的證明，且以數(shù)值算例驗(yàn)證了理論分析結(jié)果。Guo J.等[14]用L1-2公式離散Caputo導(dǎo)數(shù)項(xiàng)，用二階有限差分法離散積分項(xiàng)，得知其空間收斂階和時(shí)間收斂階都達(dá)到二階，并證明了格式的穩(wěn)定性與收斂性。Shen J.Y.等[15]提出H2N2數(shù)值微分公式（二次Hermite和Newton插值多項(xiàng)式的應(yīng)用）來(lái)近似Caputo導(dǎo)數(shù)項(xiàng)，并建立了有限差分格式。為了增加計(jì)算效率，其利用指數(shù)和近似核t1-γ，并推導(dǎo)出一種快速差分格式。對(duì)于初始時(shí)間的弱奇異性也在等級(jí)網(wǎng)格中進(jìn)行了討論。Jiang S.D.等[16]提出了分?jǐn)?shù)擴(kuò)散方程的快速L1差分格式，(0, 1)階的Caputo導(dǎo)數(shù)項(xiàng)被快速L1遞歸公式離散。Xu W.Y.等[17]考慮一種具有二階空間和時(shí)間精度的快速差分格式。用FL2-1σ公式近似時(shí)間卡普托導(dǎo)數(shù)，該公式使用了核指數(shù)和近似Caputo微分中出現(xiàn)的函數(shù)；通過(guò)離散能量方法證明了其無(wú)條件收斂性，最后通過(guò)數(shù)值例子驗(yàn)證了該方案的數(shù)值精度和效率。

本文討論如下四階分?jǐn)?shù)波方程的初邊值問(wèn)題：

其初始條件和邊界條件分別如下：

式（1）～（3）中：Ω=(0,L)×(0,T]；f(x,t)、u0(x)、ψ(x)為光滑函數(shù)；

本文中，C為一個(gè)正常數(shù)，在不同情況下可能具有不同的值。

2 預(yù)備知識(shí)

對(duì)區(qū)間[0,L]作M等分，區(qū)間[0,T]作N等分，記h=L/M,τ=T/N,xj=jh, 0≤j≤M,tn=nτ, 0≤n≤N，其中h為空間步長(zhǎng)，τ為時(shí)間步長(zhǎng)。

引理1[16]對(duì)于給定的α∈(0, 1)、截止時(shí)間限制δ、誤差ε和最后時(shí)間T，有一個(gè)正整數(shù)Nexp，正點(diǎn)和相應(yīng)的正權(quán)重，可得，其中

接下來(lái)介紹一種運(yùn)用H2N2方法計(jì)算Caputo分?jǐn)?shù)導(dǎo)數(shù)的快速算法。

令δ=τ/2，可得

其中，

引理2[15]設(shè)，γ∈(1, 2)，則有，。

其中

引理3[17]若函數(shù)，則有

3 數(shù)值離散格式

定解問(wèn)題（1）～（3）可寫成如下形式：

定義網(wǎng)格函數(shù)

將離散緊算子A分別作用于式（6）和（7），由泰勒展開式、引理1及2，可得

式（8）（9）中：

注意：初邊值條件為

引理4由帶積分余項(xiàng)的泰勒展開式及不等式（10），可得

4 緊差分格式分析

為了證明格式的收斂性，給出以下引理。

引理5[18]設(shè)u,v∈Vh，則

引理6[19]設(shè)u,v∈Vh，則有

引理7[20]設(shè)u∈Vh，可得

引理8[21]對(duì)于任意正整數(shù)p和函數(shù)G={G1,G2,…}，有

其中，{aj}定義見式（4）。

引理9[19]對(duì)任意網(wǎng)格函數(shù)，有

引理10[22]設(shè)F(k)、g(k)為非負(fù)函數(shù)，且滿足

定理1設(shè)是問(wèn)題（1）～（3）的解，是差分格式的解。為精確解，為數(shù)值解。則有

證明：令將式（8）（9）和（12）與（13）相減，得到如下誤差方程：

將算子A作用在式（13）中第二個(gè)式子，可得

利用引理2和引理3，則有

用式（16）減去式（15），可得

由引理6可知

將式（18）和（19）相加，可得

對(duì)于式（21）右端的第一項(xiàng)，利用引理9和young不等式，且考慮η0=0，可得

式（21）右端第二項(xiàng)，利用young不等式，得到

將式（22）～（25）代入式（21），有

令

因而式（26）可寫成

考慮式（10）和引理4可得

當(dāng)m為任意值時(shí)，有

由式（14）和（17）知

引用引理7，有

定理1證畢。

5 數(shù)值算例

應(yīng)用快速緊差分格式計(jì)算下列定解問(wèn)題：

其中精確解為

右端項(xiàng)為

不同步長(zhǎng)時(shí)的最大誤差為

空間收斂階為

時(shí)間收斂階為

固定時(shí)間步長(zhǎng)N=2 048，參數(shù)ε=10-12，γ取不同值時(shí)，差分格式的最大誤差以及相應(yīng)的空間收斂階與CPU時(shí)間見表1。由表1中數(shù)據(jù)可知，上述問(wèn)題的空間收斂階為4。

表1 N=2 048時(shí)最大誤差及相應(yīng)的空間階、CPU時(shí)間Table 1 Maximum error, spatial order and CPU time with N=2 048

固定空間步數(shù)M=512，參數(shù)ε=10-12，表2中列出了不同γ取值和時(shí)間步長(zhǎng)N下，計(jì)算得到的最大誤差，及其相應(yīng)的時(shí)間收斂階與CPU時(shí)間。分析表2中的數(shù)據(jù)可以得知，時(shí)間收斂階為(3-γ)階，且所需CPU時(shí)間較少。

表2 M=512時(shí)最大誤差及時(shí)間階與CPU時(shí)間Table 2 Maximum error, time order and CPU time with M=512

固定γ值，圖1給出了在空間方向上的收斂階，圖2給出了在時(shí)間方向上的收斂階。由圖1和2可看出，實(shí)例的空間收斂階與時(shí)間收斂階與理論分析得到的結(jié)論是一致的。

圖1 當(dāng)γ固定時(shí)的空間收斂階Fig.1 Spatial convergence orders with a fixed value of γ

圖2 當(dāng)γ固定時(shí)的時(shí)間收斂階Fig.2 Temporal convergence orders with a fixed value of γ

6 結(jié)語(yǔ)

本文研究了四階分?jǐn)?shù)波動(dòng)方程的快速緊差分格式，Caputo導(dǎo)數(shù)項(xiàng)用一種快速的H2N2方法來(lái)近似，通過(guò)使用降階法和離散能量法得到格式的收斂性。由數(shù)值算例可知，本文考慮的緊差分格式的收斂階為O(τ3-γ+ε+h4)，且數(shù)值算例驗(yàn)證了理論分析結(jié)果。該格式所需的CPU時(shí)間較短。