[發明專利]一種反常擴散中的基于分步計算的離散分數階差分方法

在審

• 申請號: 201710037251.0
• 申請日: 2017-01-19
• 公開（公告）號: CN106776478A
• 公開（公告）日: 2017-05-31
• 發明（設計）人: 孫洪廣;劉肖廷;達玉樂
• 申請（專利權）人: 河海大學
• 主分類號: G06F17/15
• 分類號: G06F17/15
• 代理公司: 南京縱橫知識產權代理有限公司32224
• 代理人: 朱妃,董建林

1.一種反常擴散中的基于分步計算的離散分數階差分方法，其特征在于，包括以下步驟：

1)在所考察物質的內部和邊界配置若干測試點，獲得這些測試點的擴散濃度u，確定相應的參數，包括初邊條件和擴散系數，并根據經驗公式確定分數階階數；

2)采用離散分數階差分方法離散控制方程，根據所要數值模擬物質所分布的尺度確定網格點數，網格點數包括時間點數和空間點數；其中，時間點數為N_t+1，空間點數為N_x+1，N_t為時間段數，N_x為空間段數；

3)確定反常擴散的時間空間分數階反常擴散方程，并得到時間空間反常擴散方程的離散格式；

將時間空間反常擴散方程的離散格式中的兩個gamma函數的比值看做一個參數G，分別根據時間分數階導數的階數、時間點數和空間分數階導數的階數、空間點數，計算得到時間方向gamma函數的比值G_t和空間項gamma函數的比值G_x；

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

5)通過離散格式中顯格式的數值格式，得到擴散濃度u的數值結果。

2.根據權利要求1所述的一種反常擴散中的基于分步計算的離散分數階差分方法，其特征在于：所述步驟3)中的確定反常擴散的時間空間分數階反常擴散方程，并得到時間空間反常擴散方程的離散格式，具體為，

對于描述的反常擴散的時間空間分數階反常擴散方程為，

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

時間離散分數階導數的離散形式為，

${}_0{}^C\mathrm{\Δ}{}_t{}^{\mathrm{\α}}u(x,t)=\frac{h_t^{-\mathrm{\α}}}{\mathrm{\Γ}(1-\mathrm{\α})}\underset{k=0}{\overset{N_t}{\Σ}}\frac{\mathrm{\Γ}(N_t-k+1-\mathrm{\α})}{\mathrm{\Γ}(N_t-k+1)}(u_i^{k+1}-u_i^k)---\left(2\right)$

其中，和分別為第i個點的位置第k個時刻和第k+1個時刻的濃度；擴散濃度u的上標表示時間上的時刻、下標表示空間點的位置；

空間離散分數階導數的離散形式包括右空間差分和左空間差分，

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

其中，u_n+j、u_n+j+1和u_n+j-1分別為第n+j個點、第n+j+1個點和第n+j-1個點的位置的濃度；

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

其中，u_j、u_j+1和u_j+2分別為第j個點、第j+1個點和第j+2個點的位置的濃度；

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

$u_i^{j+1}=\left(\begin{array}{c}\[h_t^{\mathrm{\α}}(D\frac{{h_x}^{-\mathrm{\β}}}{\mathrm{\Γ}(2-\mathrm{\β})}\left\{\begin{array}{c}\underset{k=0}{\overset{i}{\Σ}}\frac{\mathrm{\Γ}(k+2-\mathrm{\β})}{\mathrm{\Γ}(k+1)}(u_{i-k+1}^j-2u_{i-k}^j+u_{i-k-1}^j)\\ +\underset{k=0}{\overset{N_x-n}{\Σ}}\frac{\mathrm{\Γ}(k+2-\mathrm{\β})}{\mathrm{\Γ}(k+1)}(u_{i+k+1}^j-2u_{i+k}^j+u_{i+k-1}^j)\end{array}\right\}+A\frac{u_{i+1}^j-u_i^j}{h_x})\]\\ +\underset{k=2}{\overset{j}{\Σ}}\frac{1}{\mathrm{\Γ}(1-\mathrm{\α})}\left(\frac{\mathrm{\Γ}(j-k+1-\mathrm{\α})}{\mathrm{\Γ}(j+1-k)}-\frac{\mathrm{\Γ}(j-k+2-\mathrm{\α})}{\mathrm{\Γ}(j-k+2)}\right)u_i^k+\frac{1}{\mathrm{\Γ}(1-\mathrm{\α})}\frac{\mathrm{\Γ}(j-\mathrm{\α})}{\mathrm{\Γ}\left(j\right)}{u}_i^1\end{array}\right)---\left(5\right)$

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