本發(fā)明涉及一種離散分?jǐn)?shù)階差分方法，特別是涉及一種反常擴(kuò)散中的基于分步計(jì)算的離散分?jǐn)?shù)階差分方法，屬于數(shù)值模擬技術(shù)領(lǐng)域。

背景技術(shù)：

在環(huán)境、水利、化工、能源等許多工程領(lǐng)域，流體的動(dòng)量、熱量和質(zhì)量等基本參量在流動(dòng)過(guò)程中的運(yùn)移擴(kuò)散問(wèn)題備受關(guān)注。各類(lèi)流體流動(dòng)和傳熱傳質(zhì)問(wèn)題的根本解決，都必須基于對(duì)流體環(huán)境中運(yùn)移擴(kuò)散規(guī)律的科學(xué)把握，以便了解流體或其中含有物質(zhì)隨時(shí)空的變化情況。而隨著人類(lèi)認(rèn)知世界水平的提高，反常擴(kuò)散這個(gè)現(xiàn)象越來(lái)越普遍和受到關(guān)注。

分?jǐn)?shù)階導(dǎo)數(shù)建模已經(jīng)普遍應(yīng)用在復(fù)雜介質(zhì)中物質(zhì)的擴(kuò)散問(wèn)題，但是由于其定義的連續(xù)性，導(dǎo)致在一些非連續(xù)擴(kuò)散現(xiàn)象的建模上受到質(zhì)疑。最近，Atici，Torres等人提出了時(shí)標(biāo)上的分?jǐn)?shù)階差分，Dumitru和吳國(guó)成等人發(fā)展了離散分?jǐn)?shù)階模型，并應(yīng)用在信號(hào)處理，圖像加密和反常擴(kuò)散領(lǐng)域。但是，現(xiàn)有的離散分?jǐn)?shù)階差分方法布點(diǎn)數(shù)仍存在受限于gamma函數(shù)和計(jì)算效率不高的缺點(diǎn)。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的主要目的在于，克服現(xiàn)有技術(shù)中的不足，提供一種反常擴(kuò)散中的基于分步計(jì)算的離散分?jǐn)?shù)階差分方法，不僅擴(kuò)展原來(lái)方法模擬的范圍，減小gamma函數(shù)計(jì)算的限制，而且大幅提高計(jì)算效率，具有積極的工程實(shí)踐意義。

為了達(dá)到上述目的，本發(fā)明所采用的技術(shù)方案是：

一種反常擴(kuò)散中的基于分步計(jì)算的離散分?jǐn)?shù)階差分方法，包括以下步驟：

1)在所考察物質(zhì)的內(nèi)部和邊界配置若干測(cè)試點(diǎn)，獲得這些測(cè)試點(diǎn)的擴(kuò)散濃度u，確定相應(yīng)的參數(shù)，包括初邊條件和擴(kuò)散系數(shù)，并根據(jù)經(jīng)驗(yàn)公式確定分?jǐn)?shù)階階數(shù)；

2)采用離散分?jǐn)?shù)階差分方法離散控制方程，根據(jù)所要數(shù)值模擬物質(zhì)所分布的尺度確定網(wǎng)格點(diǎn)數(shù)，網(wǎng)格點(diǎn)數(shù)包括時(shí)間點(diǎn)數(shù)和空間點(diǎn)數(shù)；其中，時(shí)間點(diǎn)數(shù)為N_t+1，空間點(diǎn)數(shù)為N_x+1，N_t為時(shí)間段數(shù)，N_x為空間段數(shù)；

3)確定反常擴(kuò)散的時(shí)間空間分?jǐn)?shù)階反常擴(kuò)散方程，并得到時(shí)間空間反常擴(kuò)散方程的離散格式；

將時(shí)間空間反常擴(kuò)散方程的離散格式中的兩個(gè)gamma函數(shù)的比值看做一個(gè)參數(shù)G，分別根據(jù)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)、時(shí)間點(diǎn)數(shù)和空間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)、空間點(diǎn)數(shù)，計(jì)算得到時(shí)間方向gamma函數(shù)的比值G_t和空間項(xiàng)gamma函數(shù)的比值G_x；

4)將時(shí)間方向gamma函數(shù)的比值G_t和空間項(xiàng)gamma函數(shù)的比值G_x代入離散控制方程，加入初邊條件；

5)通過(guò)離散格式中顯格式的數(shù)值格式，得到擴(kuò)散濃度u的數(shù)值結(jié)果。

本發(fā)明進(jìn)一步設(shè)置為：所述步驟3)中的確定反常擴(kuò)散的時(shí)間空間分?jǐn)?shù)階反常擴(kuò)散方程，并得到時(shí)間空間反常擴(kuò)散方程的離散格式，具體為，

對(duì)于描述的反常擴(kuò)散的時(shí)間空間分?jǐn)?shù)階反常擴(kuò)散方程為，

其中，為Caputo定義下的時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)，C為Caputo定義，0和t為時(shí)間的起點(diǎn)和終點(diǎn)，α為時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)；和合在一起為Riesz定義下的空間分?jǐn)?shù)階導(dǎo)數(shù)，x^-為x位置的左鄰域，b為空間上的右端點(diǎn)，x⁺為x位置的右鄰域，a為空間上的左端點(diǎn)，β為空間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)，x^-＝(n-1)h_x,x⁺＝(n+1)h_x,x＝nh_x,1≤n≤N_x-1；x為空間上的位置，D為擴(kuò)散系數(shù)，A為對(duì)流系數(shù)，h_t為時(shí)間步長(zhǎng)，h_x為空間步長(zhǎng)，m為時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)α向上取整，v為時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)α；

時(shí)間離散分?jǐn)?shù)階導(dǎo)數(shù)的離散形式為，

其中，和分別為第i個(gè)點(diǎn)的位置第k個(gè)時(shí)刻和第k+1個(gè)時(shí)刻的濃度；擴(kuò)散濃度u的上標(biāo)表示時(shí)間上的時(shí)刻、下標(biāo)表示空間點(diǎn)的位置；

空間離散分?jǐn)?shù)階導(dǎo)數(shù)的離散形式包括右空間差分和左空間差分，

右空間差分，即右空間的離散形式為，

其中，u_n+j、u_n+j+1和u_n+j-1分別為第n+j個(gè)點(diǎn)、第n+j+1個(gè)點(diǎn)和第n+j-1個(gè)點(diǎn)的位置的濃度；

左空間差分，即左空間的離散形式為，

其中，u_j、u_j+1和u_j+2分別為第j個(gè)點(diǎn)、第j+1個(gè)點(diǎn)和第j+2個(gè)點(diǎn)的位置的濃度；

將公式(2)、公式(3)和公式(4)代入公式(1)，得到時(shí)間空間反常擴(kuò)散方程的離散格式，

考慮公式(5)中的存取在編程中很難實(shí)現(xiàn)，將調(diào)整0時(shí)刻的左邊界為其他位置的上標(biāo)下標(biāo)依次遞增。

本發(fā)明進(jìn)一步設(shè)置為：所述步驟3)中的將時(shí)間空間反常擴(kuò)散方程的離散格式中的兩個(gè)gamma函數(shù)的比值看做一個(gè)參數(shù)G，具體包括j＝1、j＝2和j＞2三種時(shí)間分?jǐn)?shù)階離散情況，以及k＝0空間分?jǐn)?shù)階離散情況。

本發(fā)明進(jìn)一步設(shè)置為：當(dāng)j＞2時(shí)，利用gamma函數(shù)的遞推公式Γ(z+1)＝zΓ(z)，結(jié)合數(shù)學(xué)中整體考慮的思想，以及時(shí)間分?jǐn)?shù)階離散結(jié)果部分，

令

其中，b＝j(luò)-k；

則

且當(dāng)j＞2時(shí)，由公式(5)可以推出b的取值范圍為(1,N_t-1)；

滿足gamma函數(shù)的遞推公式下，循環(huán)計(jì)算b從1到N_t-1的G_t的值。

本發(fā)明進(jìn)一步設(shè)置為：當(dāng)k＝0時(shí)，利用gamma函數(shù)的遞推公式Γ(z+1)＝zΓ(z)，結(jié)合數(shù)學(xué)中整體考慮的思想，以及空間分?jǐn)?shù)階離散結(jié)果部分，

令

則

由公式(5)可以推出k的取值范圍為(1,N_x-1)；

滿足gamma函數(shù)的遞推公式下，循環(huán)計(jì)算k從1到N_x-1的G_x的值。

本發(fā)明進(jìn)一步設(shè)置為：所述步驟5)中的擴(kuò)散濃度u的數(shù)值結(jié)果為，

當(dāng)j＝1時(shí)，擴(kuò)散濃度u的數(shù)值結(jié)果為，

當(dāng)j＝2時(shí)，擴(kuò)散濃度u的數(shù)值結(jié)果為，

當(dāng)j＞2時(shí)，擴(kuò)散濃度u的數(shù)值結(jié)果為，

與現(xiàn)有技術(shù)相比，本發(fā)明具有的有益效果是：

通過(guò)利用了gamma函數(shù)的遞推關(guān)系，結(jié)合數(shù)值模擬中的分布計(jì)算的思路和整體考慮的思想，將兩個(gè)gamma函數(shù)比值的問(wèn)題轉(zhuǎn)化為多個(gè)小數(shù)乘積的形式，后代入離散控制方程，得到擴(kuò)散濃度的數(shù)值結(jié)果，從而達(dá)到減小gamma函數(shù)計(jì)算的限制，擴(kuò)展模擬的點(diǎn)數(shù)，提高模擬的效率的目的，具有積極的工程實(shí)踐意義。

上述內(nèi)容僅是本發(fā)明技術(shù)方案的概述，為了更清楚的了解本發(fā)明的技術(shù)手段，下面結(jié)合附圖對(duì)本發(fā)明作進(jìn)一步的描述。

附圖說(shuō)明

圖1為本發(fā)明一種反常擴(kuò)散中的基于分步計(jì)算的離散分?jǐn)?shù)階差分方法的流程圖；

圖2為隱格式下，時(shí)刻t＝30時(shí)，不同分?jǐn)?shù)階導(dǎo)數(shù)階數(shù)在空間時(shí)間點(diǎn)數(shù)都為250情況下擴(kuò)散濃度分布情況。

具體實(shí)施方式

下面結(jié)合說(shuō)明書(shū)附圖，對(duì)本發(fā)明作進(jìn)一步的說(shuō)明。

本發(fā)明提供一種反常擴(kuò)散中的基于分步計(jì)算的離散分?jǐn)?shù)階差分方法，如圖1所示，包括以下步驟：

1)在所考察物質(zhì)的內(nèi)部和邊界配置若干測(cè)試點(diǎn)，獲得這些測(cè)試點(diǎn)的擴(kuò)散濃度u，確定相應(yīng)的參數(shù)，包括初邊條件和擴(kuò)散系數(shù)，并根據(jù)經(jīng)驗(yàn)公式確定分?jǐn)?shù)階階數(shù)。

2)采用離散分?jǐn)?shù)階差分方法離散控制方程，根據(jù)所要數(shù)值模擬物質(zhì)所分布的尺度確定網(wǎng)格點(diǎn)數(shù)，網(wǎng)格點(diǎn)數(shù)包括時(shí)間點(diǎn)數(shù)和空間點(diǎn)數(shù)；其中，時(shí)間點(diǎn)數(shù)為N_t+1，空間點(diǎn)數(shù)為N_x+1，N_t為時(shí)間段數(shù)，N_x為空間段數(shù)。

3)確定反常擴(kuò)散的時(shí)間空間分?jǐn)?shù)階反常擴(kuò)散方程，并得到時(shí)間空間反常擴(kuò)散方程的離散格式；

將時(shí)間空間反常擴(kuò)散方程的離散格式中的兩個(gè)gamma函數(shù)的比值看做一個(gè)參數(shù)G，分別根據(jù)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)、時(shí)間點(diǎn)數(shù)和空間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)、空間點(diǎn)數(shù)，計(jì)算得到時(shí)間方向gamma函數(shù)的比值G_t和空間項(xiàng)gamma函數(shù)的比值G_x；

對(duì)于描述的反常擴(kuò)散的時(shí)間空間分?jǐn)?shù)階反常擴(kuò)散方程為，

其中，為Caputo定義下的時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)，C為Caputo定義，0和t為時(shí)間的起點(diǎn)和終點(diǎn)，α為時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)；和合在一起為Riesz定義下的空間分?jǐn)?shù)階導(dǎo)數(shù)，x^-為x位置的左鄰域，b為空間上的右端點(diǎn)，x⁺為x位置的右鄰域，a為空間上的左端點(diǎn)，β為空間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)，x^-＝(n-1)h_x,x⁺＝(n+1)h_x,x＝nh_x,1≤n≤N_x-1；x為空間上的位置，D為擴(kuò)散系數(shù)，A為對(duì)流系數(shù)，h_t為時(shí)間步長(zhǎng)，h_x為空間步長(zhǎng)，m為時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)α向上取整，v為時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)α；

時(shí)間離散分?jǐn)?shù)階導(dǎo)數(shù)的離散形式為，

其中，和分別為第i個(gè)點(diǎn)的位置第k個(gè)時(shí)刻和第k+1個(gè)時(shí)刻的濃度；擴(kuò)散濃度u的上標(biāo)表示時(shí)間上的時(shí)刻、下標(biāo)表示空間點(diǎn)的位置；

空間離散分?jǐn)?shù)階導(dǎo)數(shù)的離散形式包括右空間差分和左空間差分，

右空間差分，即右空間的離散形式為，

其中，u_n+j、u_n+j+1和u_n+j-1分別為第n+j個(gè)點(diǎn)、第n+j+1個(gè)點(diǎn)和第n+j-1個(gè)點(diǎn)的位置的濃度；

左空間差分，即左空間的離散形式為，

其中，u_j、u_j+1和u_j+2分別為第j個(gè)點(diǎn)、第j+1個(gè)點(diǎn)和第j+2個(gè)點(diǎn)的位置的濃度；

將公式(2)、公式(3)和公式(4)代入公式(1)，得到時(shí)間空間反常擴(kuò)散方程的離散格式，

考慮公式(5)中的存取在編程中很難實(shí)現(xiàn)，將調(diào)整0時(shí)刻的左邊界為其他位置的上標(biāo)下標(biāo)依次遞增。考慮公式(5)中的存取在編程中很難實(shí)現(xiàn)，將調(diào)整0時(shí)刻的左邊界為其他位置的上標(biāo)下標(biāo)依次遞增。

將時(shí)間空間反常擴(kuò)散方程的離散格式中的兩個(gè)gamma函數(shù)的比值看做一個(gè)參數(shù)G，具體包括j＝1、j＝2和j＞2三種時(shí)間分?jǐn)?shù)階離散情況，以及k＝0空間分?jǐn)?shù)階離散情況。

當(dāng)j＞2時(shí)，利用gamma函數(shù)的遞推公式Γ(z+1)＝zΓ(z)，結(jié)合數(shù)學(xué)中整體考慮的思想，以及時(shí)間分?jǐn)?shù)階離散結(jié)果部分，

令

其中，b＝j(luò)-k；

則

且當(dāng)j＞2時(shí)，由公式(5)可以推出b的取值范圍為(1,N_t-1)；

滿足gamma函數(shù)的遞推公式下，循環(huán)計(jì)算b從1到N_t-1的G_t的值。

當(dāng)k＝0時(shí)，利用gamma函數(shù)的遞推公式Γ(z+1)＝zΓ(z)，結(jié)合數(shù)學(xué)中整體考慮的思想，以及空間分?jǐn)?shù)階離散結(jié)果部分，

令

則

由公式(5)可以推出k的取值范圍為(1,N_x-1)；

滿足gamma函數(shù)的遞推公式下，循環(huán)計(jì)算k從1到N_x-1的G_x的值。

4)將時(shí)間方向gamma函數(shù)的比值G_t和空間項(xiàng)gamma函數(shù)的比值G_x代入離散控制方程，加入初邊條件；

5)通過(guò)離散格式中顯格式的數(shù)值格式，得到擴(kuò)散濃度u的數(shù)值結(jié)果。

當(dāng)j＝1時(shí)，擴(kuò)散濃度u的數(shù)值結(jié)果為，

當(dāng)j＝2時(shí)，擴(kuò)散濃度u的數(shù)值結(jié)果為，

當(dāng)j＞2時(shí)，擴(kuò)散濃度u的數(shù)值結(jié)果為，

實(shí)施例：

考慮如下離散分?jǐn)?shù)階泥沙輸運(yùn)問(wèn)題，控制方程為：

其中，擴(kuò)散系數(shù)D＝0.04，對(duì)流系數(shù)A＝0.05，時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)α＝0.8，空間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)β＝1.6。

為了驗(yàn)證方法的有效性，我們選取不同時(shí)間步長(zhǎng)和空間步長(zhǎng)來(lái)驗(yàn)證方法是否減小了gamma函數(shù)的限制和提高計(jì)算效率，對(duì)流項(xiàng)采用隱格式有限差分算法離散。

精確數(shù)值結(jié)果為如下：

對(duì)于j＝1的情況，擴(kuò)散濃度u的數(shù)值結(jié)果為，

對(duì)于j＝2的情況，擴(kuò)散濃度u的數(shù)值結(jié)果為，

對(duì)于j＞2的情況，擴(kuò)散濃度u的數(shù)值結(jié)果為，

在運(yùn)行條件為Inter(R)Core(TM)i5-5490CPU@3.30GHz下，并在隱格式下，統(tǒng)計(jì)了控制方程不同時(shí)間和空間點(diǎn)數(shù)的計(jì)算時(shí)間，無(wú)法計(jì)算則沒(méi)有結(jié)果，如表1所示。

從表1中可以發(fā)現(xiàn)，本發(fā)明提供的一種反常擴(kuò)散中的基于分步計(jì)算的離散分?jǐn)?shù)階差分方法，其計(jì)算范圍減少了gamma函數(shù)的限制，而且和已有的離散分?jǐn)?shù)階差分方法相比，該方法可以提高一倍的計(jì)算效率。

表1

圖2顯示了該方法隱格式下，在不同的時(shí)間和空間分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)條件下的時(shí)間空間布點(diǎn)數(shù)為250的計(jì)算結(jié)果，顯示了該方法在反常擴(kuò)散模擬的結(jié)果。從圖2可以得出，本發(fā)明所提供的一種反常擴(kuò)散中的基于分步計(jì)算的離散分?jǐn)?shù)階差分方法具有好的效率和擴(kuò)展離散分?jǐn)?shù)階模型應(yīng)用范圍的效果。

本發(fā)明的創(chuàng)新點(diǎn)在于，不增加計(jì)算量和存儲(chǔ)量，計(jì)算簡(jiǎn)單；將兩個(gè)gamma函數(shù)比值的問(wèn)題轉(zhuǎn)化為多個(gè)小數(shù)乘積的形式，后代入離散的擴(kuò)散方程，從而達(dá)到減小gamma函數(shù)計(jì)算的限制，擴(kuò)展了模擬的點(diǎn)數(shù)，提高了模擬的效率；該方法簡(jiǎn)單快速，易于實(shí)施，方便工程技術(shù)人員的使用，可用模擬污染物在土壤中的遷移、泥沙輸運(yùn)、石油滲流、地下水的傳輸和湍流等水分和溶質(zhì)遷移擴(kuò)散過(guò)程，為科學(xué)的防止和治理環(huán)境污染提供了科學(xué)的指導(dǎo)，具有重要的理論和工程意義。

以上顯示和描述了本發(fā)明的基本原理、主要特征及優(yōu)點(diǎn)。本行業(yè)的技術(shù)人員應(yīng)該了解，本發(fā)明不受上述實(shí)施例的限制，上述實(shí)施例和說(shuō)明書(shū)中描述的只是說(shuō)明本發(fā)明的原理，在不脫離本發(fā)明精神和范圍的前提下，本發(fā)明還會(huì)有各種變化和改進(jìn)，這些變化和改進(jìn)都落入要求保護(hù)的本發(fā)明范圍內(nèi)。本發(fā)明要求保護(hù)范圍由所附的權(quán)利要求書(shū)及其等效物界定。