matlab和Python常微分方程的非標準有限差分法
標準有限差分方案
遇到標準有限差分法的主要問題是數值不穩定性。如果存在有限差分方程的解與微分方程的任何可能解都不對應，則稱微分方程的離散模型具有數值不穩定性。 Mickens指出了一些數值不穩定原因。這些原因是：
-
用不同階數的離散導數表示微分方程中的導數。 這種情況的一個例子是在表示邏輯模型時
d P d t = P ( 1 − P ) , P ( 0 ) = P 0 \frac{\boldsymbol{dP}}{\boldsymbol{dt}}=\boldsymbol{P}\left( 1-\boldsymbol{P} \right) ,\boldsymbol{P}\left( 0 \right) =\boldsymbol{P}_0
dtdP=P(1−P),P(0)=P0使用中心差分方案：
P j + 1 − P j − 1 2 h = P j ( 1 − P j ) ⇒ P j + 1 = P j − 1 + 2 h P j ( 1 − P j ) \frac{\boldsymbol{P}^{\boldsymbol{j}+1}-\boldsymbol{P}^{\boldsymbol{j}-1}}{2\boldsymbol{h}}=\boldsymbol{P}^{\boldsymbol{j}}\left( 1-\boldsymbol{P}^{\boldsymbol{j}} \right) \Rightarrow \boldsymbol{P}^{\boldsymbol{j}+1}=\boldsymbol{P}^{\boldsymbol{j}-1}+2\boldsymbol{hP}^{\boldsymbol{j}}\left( 1-\boldsymbol{P}^{\boldsymbol{j}} \right)
2hPj+1−Pj−1=Pj(1−Pj)⇒Pj+1=Pj−1+2hPj(1−Pj)下圖1顯示了使用中心，前向和後向差分方案的邏輯模型 P ˙ ( t ) = P ( t ) ( 1 − P ( t ) ) , P ( 0 ) = 0.5 \boldsymbol{\dot{P}}\left( \boldsymbol{t} \right) =\boldsymbol{P}\left( \boldsymbol{t} \right) \left( 1-\boldsymbol{P}\left( \boldsymbol{t} \right) \right) ,\boldsymbol{P}\left( 0 \right) =0.5
P˙(t)=P(t)(1−P(t)),P(0)=0.5的解。 該圖顯示，通過前向和後向差分方案獲得的解與原始模型具有相同的行為，而中心有限差分方案與精確解無關。圖1
-
使用標準分母函式 h \boldsymbol{h} h。 在點 t j \boldsymbol{t}_{\boldsymbol{j}} tj的一階和二階導數近似為：
d u d t ( t i ) ≈ u ( t i + 1 ) − u ( t i ) h 或 d u d t ( t i ) ≈ u ( t i ) − u ( t i − 1 ) h d 2 u d t 2 ( t i ) ≈ u ( t i + 1 ) − 2 u ( t i ) + u ( t i + 1 ) h 2 \frac{\boldsymbol{du}}{\boldsymbol{dt}}\left( \boldsymbol{t}_{\boldsymbol{i}} \right) \approx \frac{\boldsymbol{u}\left( \boldsymbol{t}_{\boldsymbol{i}+1} \right) -\boldsymbol{u}\left( \boldsymbol{t}_{\boldsymbol{i}} \right)}{\boldsymbol{h}}或\frac{\boldsymbol{du}}{\boldsymbol{dt}}\left( \boldsymbol{t}_{\boldsymbol{i}} \right) \approx \frac{\boldsymbol{u}\left( \boldsymbol{t}_{\boldsymbol{i}} \right) -\boldsymbol{u}\left( \boldsymbol{t}_{\boldsymbol{i}-1} \right)}{\boldsymbol{h}}\\\frac{\boldsymbol{d}^2\boldsymbol{u}}{\boldsymbol{dt}^2}\left( \boldsymbol{t}_{\boldsymbol{i}} \right) \approx \frac{\boldsymbol{u}\left( \boldsymbol{t}_{\boldsymbol{i}+1} \right) -2\boldsymbol{u}\left( \boldsymbol{t}_{\boldsymbol{i}} \right) +\boldsymbol{u}\left( \boldsymbol{t}_{\boldsymbol{i}+1} \right)}{\boldsymbol{h}^2} dtdu(ti)≈hu(ti+1)−u(ti)或dtdu(ti)≈hu(ti)−u(ti−1)dt2d2u(ti)≈h2u(ti+1)−2u(ti)+u(ti+1)
在許多情況下，經典分母函式的選擇不合適，並且可能導致數值不穩定。 請考慮下面的示例，一階IVP：
d y ( t ) d t = A y ( t ) + t , y ( 0 ) = 0 \frac{\boldsymbol{dy}\left( \boldsymbol{t} \right)}{\boldsymbol{dt}}=\boldsymbol{Ay}\left( \boldsymbol{t} \right) +\boldsymbol{t},\boldsymbol{y}\left( 0 \right) =0 dtdy(t)=Ay(t)+t,y(0)=0
這個IVP確切解是
y ( t ) = 1 A 2 ( e A t − ( A t + 1 ) ) \boldsymbol{y}\left( \boldsymbol{t} \right) =\frac{1}{\boldsymbol{A}^2}\left( \boldsymbol{e}^{\boldsymbol{At}}-\left( \boldsymbol{At}+1 \right) \right) y(t)=A21(eAt−(At+1))
間隔 [ 0 , 1 ] \left[ 0,1 \right] [0,1] 將被劃分為1000個子間隔，按照點 0 = t 0 < t 1 < ⋯ < t 1000 = 1 0=\boldsymbol{t}_0<\boldsymbol{t}_1<\cdots <\boldsymbol{t}_{1000}=1 0=t0<t1<⋯<t1000=1，每個子間隔的長度 h = 0.001 \boldsymbol{h}=0.001 h=0.001。 前向尤拉方法，隱式梯形法則和經典四階龍格庫塔法將用於在離散點處針對引數A的不同值解決給定問題。從理論上講，尤拉方法的誤差應為 O ( h ) = O ( 1 0 − 3 ) \mathcal{O}\left( \boldsymbol{h} \right) =\mathcal{O}\left( 10^{-3} \right) O(h)=O(10−3) ，隱式梯形規則的誤差應為 O ( h 2 ) = O ( 1 0 − 6 ) \mathcal{O}\left( \boldsymbol{h}^2 \right) =\mathcal{O}\left( 10^{-6} \right) O(h2)=O(10−6) ，而四階龍格庫塔法誤差應為 O ( h 4 ) = O ( 1 0 − 12 ) \mathcal{O}\left( \boldsymbol{h}^4 \right) =\mathcal{O}\left( 10^{-12} \right) O(h4)=O(10−12) 。
分別計算三種方法的不同A值誤差的Python程式碼是：
執行程式碼可獲得以下結果：
runfile(’/media/WindowsData1/PyFiles/ numinstduetotrivstepfun.py’, wdir=’/media/WindowsData1/PyFiles’) -------------------------------------------------------- A Euler Error Imp. Trapz Err RK4 Error -------------------------------------------------------- 1 1.35789622e-03 3.56321584e-07 2.30926389e-14 3 1.00002555e-02 5.19127795e-06 4.50417481e-12 5 7.35013397e-02 4.52531609e-05 1.53952406e-10 7 5.39170234e-01 3.52721966e-04 3.11633030e-09 9 3.94740396e+00 2.65343074e-03 4.88584107e-08 11 2.88444591e+01 1.97153658e-02 6.58044428e-07 13 2.10373074e+02 1.45836606e-01 8.01259330e-06 15 1.53146894e+03 1.07702617e+00 9.07992144e-05 17 1.11283670e+04 7.94936660e+00 9.75035611e-04 19 8.07191089e+04 5.86610000e+01 1.00415309e-02 21 5.84468245e+05 4.32849955e+02 1.00014394e-01 23 4.22478467e+06 3.19388027e+03 9.69289817e-01 25 3.04879507e+07 2.35668600e+04 9.18238944e+00 27 2.19661624e+08 1.73896191e+05 8.53281952e+01 29 1.58017839e+09 1.28317310e+06 7.79939406e+02 --------------------------------------------------------可以看出，儘管步長 h \boldsymbol{h} h保持恆定，但是三種方法的誤差都隨著 A \boldsymbol{A} A的增加而增加。
-
使用非線性項的區域性逼近：標準有限差分方案使用諸如 x j 2 或 x j + 1 2 \boldsymbol{x}^{\boldsymbol{j}^2}或\boldsymbol{x}^{\boldsymbol{j}+1^2} xj2或xj+12的區域性逼近來表示諸如 x 2 ( t j ) \boldsymbol{x}^2\left( \boldsymbol{t}_{\boldsymbol{j}} \right) x2(tj) 的非線性項。 Mickens表明，非線性項的這種區域性逼近可能導致數值不穩定。 Mickens給出的一個例子是指數衰減模型的解：
d P ( t ) d t = − P ( t ) , P ( 0 ) = 0.5 \frac{\boldsymbol{dP}\left( \boldsymbol{t} \right)}{\boldsymbol{dt}}=-\boldsymbol{P}\left( \boldsymbol{t} \right) ,\boldsymbol{P}\left( 0 \right) =0.5 dtdP(t)=−P(t),P(0)=0.5
指數衰減模型的前向尤拉離散化由下式給出：
P n + 1 = P n − h P n = ( 1.0 − h ) P n = ( 1 − h ) n P 0 \boldsymbol{P}^{\boldsymbol{n}+1}=\boldsymbol{P}^{\boldsymbol{n}}-\boldsymbol{hP}^{\boldsymbol{n}}=\left( 1.0-\boldsymbol{h} \right) \boldsymbol{P}^{\boldsymbol{n}}=\left( 1-\boldsymbol{h} \right) ^{\boldsymbol{n}}\boldsymbol{P}^0 Pn+1=Pn−hPn=(1.0−h)Pn=(1−h)nP0
對於步長 h h h的不同值，所得結果如下圖2所示。
圖2
為指數衰減模型獲得的解包括具有漸近穩定，週期性和不穩定動力學的解曲線，其中週期性和不穩定解不對應於指數衰減模型的任何真實解。 僅人口圖衰減為零的第一個圖就與微分方程的真實解一致。
另一個例子是使用亨恩的方法來求解邏輯模型，
d x ( t ) d t = x ( t ) ( 1 − x ( t ) ) , x ( 0 ) = 0.5 \frac{\boldsymbol{dx}\left( \boldsymbol{t} \right)}{\boldsymbol{dt}}=\boldsymbol{x}\left( \boldsymbol{t} \right) \left( 1-\boldsymbol{x}\left( \boldsymbol{t} \right) \right) ,\boldsymbol{x}\left( 0 \right) =0.5 dtdx(t)=x(t)(1−x(t)),x(0)=0.5
Python程式碼如下，
下圖3顯示了針對亨恩的方法對 h \boldsymbol{h} h的不同值進行邏輯模型的求解。
圖3
從上圖3可以看出，與步長 h \boldsymbol{h} h的某些值相對應，邏輯模型的數值解可能具有周期性和混沌行為，這與模型的任何解都不對應。
非標準有限差分格式的構造規則
與標準有限差分方案一起設計的數值穩定性表明，由於標準模型中的某些固有誤差，檢索真實資訊變得非常昂貴，導致這些標準方案無法準確顯示。
考慮到導致標準有限差分方案中數值不穩定的原因，Micken為這種數值方案的設計建立了四個非標準構造規則。 這些規則為：
-
因為使用與微分方程的階不同的離散導數會導致數值不穩定，所以離散模型的階應與微分方程的階相同。 在此規則下，中心有限差分方案不能用作一階微分方程離散模型中一階導數的近似值。 可以使用前向或後向差分方案。
-
非標準分母函式必須用於連續導數的離散表示。 d y d t \frac{\boldsymbol{dy}}{\boldsymbol{dt}} dtdy在 t = t j \boldsymbol{t}=\boldsymbol{t}_{\boldsymbol{j}} t=tj的非標準離散表示形式為：
d y d t ≈ y ( t j + 1 ) − ψ ( λ , h ) y ( t j ) ϕ ( λ , h ) \frac{\boldsymbol{dy}}{\boldsymbol{dt}}\approx \frac{\boldsymbol{y}\left( \boldsymbol{t}_{\boldsymbol{j}+1} \right) -\boldsymbol{\psi }\left( \boldsymbol{\lambda },\boldsymbol{h} \right) \boldsymbol{y}\left( \boldsymbol{t}_{\boldsymbol{j}} \right)}{\boldsymbol{\phi }\left( \boldsymbol{\lambda },\boldsymbol{h} \right)} dtdy≈ϕ(λ,h)y(tj+1)−ψ(λ,h)y(tj)
其中 λ \boldsymbol{\lambda } λ是模型引數的向量， h = t j + 1 − t j \boldsymbol{h}=\boldsymbol{t}_{\boldsymbol{j}+1}-\boldsymbol{t}_{\boldsymbol{j}} h=tj+1−tj。 分子和分母函式 ψ ( λ , h ) \boldsymbol{\psi }\left( \boldsymbol{\lambda },\boldsymbol{h} \right) ψ(λ,h) 和 ϕ ( λ , h ) \boldsymbol{\phi }\left( \boldsymbol{\lambda },\boldsymbol{h} \right) ϕ(λ,h) 應滿足以下特性：
當 h → 0 , ψ ( λ , h ) → 1 與 ϕ ( λ , h ) → h 當\boldsymbol{h}\rightarrow 0,\boldsymbol{\psi }\left( \boldsymbol{\lambda },\boldsymbol{h} \right) \rightarrow 1與\boldsymbol{\phi }\left( \boldsymbol{\lambda },\boldsymbol{h} \right) \rightarrow \boldsymbol{h} 當h→0,ψ(λ,h)→1與ϕ(λ,h)→h
舉例說明使用非平凡分母如何導致數值不穩定，請考慮指數衰減模型
d P ( t ) d t = − P ( t ) , P ( 0 ) = 0.5 \frac{\boldsymbol{dP}\left( \boldsymbol{t} \right)}{\boldsymbol{dt}}=-\boldsymbol{P}\left( \boldsymbol{t} \right) ,\boldsymbol{P}\left( 0 \right) =0.5 dtdP(t)=−P(t),P(0)=0.5
在尤拉方法中，假設平凡的分母函式 h \boldsymbol{h} h被分母函式代替
ϕ ( h ) = 1 − e − h h \boldsymbol{\phi }\left( \boldsymbol{h} \right) =\frac{1-\boldsymbol{e}^{-\boldsymbol{h}}}{\boldsymbol{h}} ϕ(h)=h1−e−h
因此，所得離散模型為
P j + 1 = ( 1 − ϕ ) P j , P 0 = 0.5 \boldsymbol{P}^{\boldsymbol{j}+1}=\left( 1-\boldsymbol{\phi } \right) \boldsymbol{P}^{\boldsymbol{j}},\boldsymbol{P}^0=0.5 Pj+1=(1−ϕ)Pj,P0=0.5
Python程式碼expdecnsden.py用於求解步長 h \boldsymbol{h} h的不同值的指數衰減模型：
執行該程式碼的結果如下圖4所示，其中對於不同步長 h \boldsymbol{h} h的值，使用非標準分母函式求解了指數衰減模型。
-
非線性項應使用非區域性近似來近似：像 y 2 ( t ) \boldsymbol{y}^2\left( \boldsymbol{t} \right) y2(t) 的項應由 y j y j + 1 \boldsymbol{y}^{\boldsymbol{j}}\boldsymbol{y}^{\boldsymbol{j}+1} yjyj+1或 y j − 1 y j + 1 \boldsymbol{y}^{\boldsymbol{j}-1}\boldsymbol{y}^{\boldsymbol{j}+1} yj−1yj+1代替 ( y j ) 2 \left( \boldsymbol{y}^{\boldsymbol{j}} \right) ^2 (yj)2來表示，其中 ( y j ) 2 \left( \boldsymbol{y}^{\boldsymbol{j}} \right) ^2 (yj)2。 例如，對應於邏輯模型的離散模型
d y d t = y ( t ) ( 1 − y ( t ) ) , y ( 0 ) = 0.5 \frac{\boldsymbol{dy}}{\boldsymbol{dt}}=\boldsymbol{y}\left( \boldsymbol{t} \right) \left( 1-\boldsymbol{y}\left( \boldsymbol{t} \right) \right) ,\boldsymbol{y}\left( 0 \right) =0.5 dtdy=y(t)(1−y(t)),y(0)=0.5
可能是
y j + 1 − y j h = y j ( 1 − y j + 1 ) ⇒ y j + 1 = y j ( 1 + h ) 1 + h y j , y 0 = 0.5 \frac{\boldsymbol{y}^{\boldsymbol{j}+1}-\boldsymbol{y}^{\boldsymbol{j}}}{\boldsymbol{h}}=\boldsymbol{y}^{\boldsymbol{j}}\left( 1-\boldsymbol{y}^{\boldsymbol{j}+1} \right) \Rightarrow \boldsymbol{y}^{\boldsymbol{j}+1}=\frac{\boldsymbol{y}^{\boldsymbol{j}}\left( 1+\boldsymbol{h} \right)}{1+\boldsymbol{hy}^{\boldsymbol{j}}},\boldsymbol{y}^0=0.5 hyj+1−yj=yj(1−yj+1)⇒yj+1=1+hyjyj(1+h),y0=0.5
Python程式碼nonlocapplog.py用於使用非區域性逼近法求解邏輯模型。 該程式碼是
下圖5是通過執行python程式碼獲得的。 它顯示了使用不同步長值的邏輯模型的解決方案。 所有這些解都對應於微分方程模型的真實解。 因此，對於步長的所有值，離散模型不會遭受數值不穩定。
圖5
上圖5是通過執行python程式碼獲得的。 它顯示了使用不同步長值的邏輯模型的解決方案。 所有這些解都對應於微分方程模型的真實解。 因此，對於步長的所有值，離散模型不會遭受數值上的不穩定性。
-
離散模型和微分方程的解之間的動態一致性。 這意味著微分方程解所具有的所有性質都必須由離散模型的解所具有。 特殊屬性包括正性，單調性，有界性，極限環和其他週期解等。
非標準有限差分方案是基於非標準規則構造的微分方程的離散表示。