隱式顯式方法避免了直接求解非線性問題。這對某些問題可能是有利的，但在其他問題中可能導致嚴重的步長限制。此外，由此產生的數值方案有時可能具有不良的定性特性。出於這個原因，我們需要描述允許我們求解完全隱式數值方案中產生的非線性方程的方法。
我們考慮一個常微分方程 d y d t = f ( t , y ) {\displaystyle {\frac {\mathrm {d} y}{\mathrm {d} t}}=f(t,y)} t ∈ [ t 0 , t ∗ ] {\displaystyle t\in [t_{0},t^{*}]} f ( t , y ) {\displaystyle f(t,y)} y {\displaystyle y}
以下結果主要來自 Iserles[ 1] [ 2] [ 3] t i {\displaystyle t_{i}} i {\displaystyle i} y i {\displaystyle y_{i}} i {\displaystyle i} h {\displaystyle h} f {\displaystyle f}
我們將證明這兩種方法的收斂性（修正牛頓-拉夫森方法的收斂性證明見 Iserles[ 4] y ′ ( t ) = y 2 {\displaystyle y'(t)=y^{2}} y ( 0 ) = 1 {\displaystyle y(0)=1} t ∈ [ 0 , 0.99 ] {\displaystyle t\in [0,0.99]}
我們考慮 d y d t = y 2 {\displaystyle {\frac {\mathrm {d} y}{\mathrm {d} t}}=y^{2}} y ( t = 0 ) = 1 {\displaystyle y(t=0)=1} t ∈ [ 0 , 0.99 ] {\displaystyle t\in [0,0.99]} y ( t ) {\displaystyle y(t)} d y d t {\displaystyle {\frac {\mathrm {d} y}{\mathrm {d} t}}}
為了顯式地積分這個方程，我們使用變數分離法，發現 ∫ y ( 0 ) y ( t ) 1 y ~ 2 d y ~ = ∫ 0 t d τ {\displaystyle {}\int _{y(0)}^{y(t)}{\frac {1}{{\tilde {y}}^{2}}}\mathrm {d} {\tilde {y}}=\int _{0}^{t}\mathrm {d} \tau } − 1 y ( t ) = t + c {\displaystyle -{\frac {1}{y(t)}}=t+c} c {\displaystyle c} c = − 1 {\displaystyle c=-1} y ( t ) = 1 1 − t {\displaystyle y(t)={\frac {1}{1-t}}} t → 1 {\displaystyle t\rightarrow 1}
定義 函式 f ( x ) : x ∈ D ⊂ R {\displaystyle f(x):x\in D\subset \mathbb {R} } 利普希茨 函式，如果對於所有 x 1 {\displaystyle x_{1}} x 2 {\displaystyle x_{2}} D {\displaystyle D} ‖ f ( x 1 ) − f ( x 2 ) ‖ ≤ λ ‖ x 1 − x 2 ‖ {\displaystyle \left\Vert f(x_{1})-f(x_{2})\right\Vert \leq \lambda \left\Vert x_{1}-x_{2}\right\Vert }
利普希茨條件有兩個具體的定義。
定義 函式 f ( x ) {\displaystyle f(x)} 區域性利普希茨 函式，如果對於每個 z ∈ R {\displaystyle z\in \mathbb {R} } L > 0 {\displaystyle L>0} f {\displaystyle f} z {\displaystyle z} L {\displaystyle L}
定義 如果 f ( x ) {\displaystyle f(x)} R {\displaystyle \mathbb {R} } R {\displaystyle \mathbb {R} } f {\displaystyle f} 全域性利普希茨 函式。
請注意區域性利普希茨條件和全域性利普希茨條件之間的根本區別。在區域性版本中，利普希茨“常數”（ λ ) {\displaystyle \lambda )} x ∈ R {\displaystyle x\in \mathbb {R} } λ ) {\displaystyle \lambda )} R {\displaystyle \mathbb {R} }
Peano 定理指出，如果 f ( x ) {\displaystyle f(x)} x ′ ( t ) = f ( x ) {\displaystyle x'(t)=f(x)} t 0 {\displaystyle t_{0}} f ( x ) {\displaystyle f(x)} x ′ ( t ) = f ( x ) {\displaystyle x'(t)=f(x)} x ( t 0 ) = x 0 {\displaystyle x(t_{0})=x_{0}} [ 5] [ 6]
我們回憶一下後向尤拉方法由 y n + 1 = y n + h f ( y n + 1 ) . {\displaystyle y^{n+1}=y^{n}+hf(y^{n+1}).} f {\displaystyle f} y n + 1 = y n + h ( y n + 1 ) 2 {\displaystyle y^{n+1}=y^{n}+h\left(y^{n+1}\right)^{2}}
我們經常使用函式迭代來求解非線性方程。我們回憶一下，有兩種常用的方法：不動點迭代和牛頓法。
我們想要找到一個方程 x = f ( x ) . {\displaystyle x=f(x).} x n + 1 = f ( x n ) {\displaystyle x_{n+1}=f(x_{n})} n = 0 , 1 , 2 … {\displaystyle n=0,1,2\ldots }
定理 設 f ( x ) {\displaystyle f(x)} x ~ = f ( x ~ ) {\displaystyle {\tilde {x}}=f({\tilde {x}})} x ∈ ( a , b ) ⊂ R {\displaystyle x\in (a,b)\subset \mathbb {R} } k < 1 {\displaystyle k<1} f ( x ) {\displaystyle f(x)} [ a , b ] {\displaystyle [a,b]} x n + 1 = f ( x n ) {\displaystyle x_{n+1}=f(x_{n})} x ~ = x ∞ = f ( x ∞ ) {\displaystyle {\tilde {x}}=x_{\infty }=f(x_{\infty })} x ∈ [ a , b ] {\displaystyle x\in [a,b]} f ( x ) {\displaystyle f(x)}
證明 由於 f ( x ) {\displaystyle f(x)} | x n + 1 − x ∞ | = | f ( x n ) − f ( x ∞ ) | ≤ k | x n − x ∞ | {\displaystyle \left|x_{n+1}-x_{\infty }\right|=\left|f(x_{n})-f(x_{\infty })\right|\leq k\left|x_{n}-x_{\infty }\right|} n = 1 , 2 … {\displaystyle n=1,2\ldots } | x n + 1 − x ∞ | ≤ k n | x 1 − x ∞ | . {\displaystyle \left|x_{n+1}-x_{\infty }\right|\leq k^{n}\left|x_{1}-x_{\infty }\right|.} k < 1 {\displaystyle k<1} lim n → ∞ k n | x 1 − x ∞ | = 0 {\displaystyle \lim _{n\rightarrow \infty }k^{n}|x_{1}-x_{\infty }|=0} x ∞ = f ( x ∞ ) {\displaystyle x_{\infty }=f(x_{\infty })} x ∞ {\displaystyle x_{\infty }} ◻ {\displaystyle \Box }
有關在該定理中使用的假設下不動點存在的證明，請參閱數值分析方面的書籍，例如 Bradie[ 7] [ 8]
關於我們的問題，我們應用不動點迭代，我們想要找到以下形式方程的根
w = h w 2 + β = f ( w ) . {\displaystyle w=hw^{2}+\beta =f(w).}
當時間步長 h {\displaystyle h} f ′ ( w ) = 2 h w ≤ 200 h < 1 {\displaystyle f'(w)=2hw\leq 200h<1} 1
現在我們考慮牛頓法。我們想要找到一個根， x ∗ {\displaystyle x^{*}} f ( x ) {\displaystyle f(x)} f ( x ∗ ) = 0. {\displaystyle f(x^{*})=0.}
x n + 1 = x n − f ( x n ) f ′ ( x n ) {\displaystyle x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}}
其中 n = 0 , 1 , 2 … {\displaystyle n=0,1,2\ldots } f ( x ) {\displaystyle f(x)}
定理 假設 f ( x ) {\displaystyle f(x)} x ∗ {\displaystyle x^{*}} f ( x ∗ ) = 0 {\displaystyle f(x^{*})=0} f ′ ( x ) ≠ 0 {\displaystyle f'(x)\neq 0} [ x ∗ − | x ∗ − x 0 | , x ∗ + | x ∗ − x 0 | ] {\displaystyle \left[x^{*}-\left|x^{*}-x_{0}\right|,x^{*}+\left|x^{*}-x_{0}\right|\right]} f ″ ( x ) {\displaystyle f''(x)} | x 0 − x ∗ | {\displaystyle \left|x_{0}-x^{*}\right|}
Proof f ( x ∗ ) = f ( x n ) + f ′ ( x n ) ( x ∗ − x n ) + 1 2 ! f ″ ( z n ) ( x ∗ − x n ) 2 {\displaystyle f(x^{*})=f(x_{n})+f'(x_{n})(x^{*}-x_{n})+{\frac {1}{2!}}f''(z_{n})(x^{*}-x_{n})^{2}} z n ∈ [ x n , x ∗ ] {\displaystyle z_{n}\in [x_{n},x^{*}]} f ( x ∗ ) = 0 {\displaystyle f(x^{*})=0} 0 = f ( x n ) + f ′ ( x n ) ( x ∗ − x n ) + 1 2 ! f ″ ( z n ) ( x ∗ − x n ) 2 , {\displaystyle 0=f(x_{n})+f'(x_{n})(x^{*}-x_{n})+{\frac {1}{2!}}f''(z_{n})(x^{*}-x_{n})^{2},} f ( x n ) f ′ ( x n ) + ( x ∗ − x n ) = − 1 2 ! f ″ ( z n ) f ′ ( z n ) ( x ∗ − x n ) 2 . {\displaystyle {\frac {f(x_{n})}{f'(x_{n})}}+(x^{*}-x_{n})=-{\frac {1}{2!}}{\frac {f''(z_{n})}{f'(z_{n})}}(x^{*}-x_{n})^{2}.} x n + 1 {\displaystyle x_{n+1}} 2 x ∗ − x n + 1 = − 1 2 ! f ″ ( z n ) f ′ ( z n ) ( x ∗ − x n ) 2 . {\displaystyle x^{*}-x_{n+1}=-{\frac {1}{2!}}{\frac {f''(z_{n})}{f'(z_{n})}}(x^{*}-x_{n})^{2}.} e n = | x ∗ − x n | . {\displaystyle e_{n}=\left|x^{*}-x_{n}\right|.} e n + 1 = | 1 2 ! f ″ ( z n ) f ′ ( z n ) | e n 2 {\displaystyle e_{n+1}=\left|{\frac {1}{2!}}{\frac {f''(z_{n})}{f'(z_{n})}}\right|e_{n}^{2}} c {\displaystyle c} | 1 2 ! f ″ ( z n ) f ′ ( z n ) | < c . {\displaystyle \left|{\frac {1}{2!}}{\frac {f''(z_{n})}{f'(z_{n})}}\right|<c.} e n + 1 < m e n 2 {\displaystyle e_{n+1}<me_{n}^{2}} m {\displaystyle m} e 0 = | x 0 − x ∗ | {\displaystyle e_{0}=|x_{0}-x^{*}|} ◻ {\displaystyle \Box }
關於我們的問題，我們考慮 f ( y ) = y − h y 2 − β {\displaystyle f(y)=y-hy^{2}-\beta } f ′ ( y ) = 1 − 2 h y ≠ 0 {\displaystyle f'(y)=1-2hy\neq 0} f ″ ( y ) {\displaystyle f''(y)}
向後尤拉法、向前尤拉法和克朗尼克森法是 Theta 方法的特殊情況，因此我們將首先證明 Theta 方法的收斂性以包含這三種方法。Theta 方法是以下演算法： y n + 1 = y n + h [ θ f ( t n , y n ) + ( 1 − θ ) f ( t n + 1 , y n + 1 ) ] {\displaystyle y^{n+1}=y^{n}+h[\theta f(t^{n},y^{n})+(1-\theta )f(t^{n+1},y^{n+1})]} n = 0 , 1 , … {\displaystyle n=0,1,\ldots } θ ∈ [ 0 , 1 ] {\displaystyle \theta \in [0,1]} θ = 1 / 2 {\displaystyle \theta =1/2}
首先，代入精確解 y ( t ) {\displaystyle y(t)}
y ( t n + 1 ) − y ( t n ) − h [ θ f ( t n , y ( t n ) ) + ( 1 − θ ) f ( t n + 1 , y ( t n + 1 ) ) ] {\displaystyle {}y(t^{n+1})-y(t^{n})-h[\theta f(t^{n},y(t^{n}))+(1-\theta )f(t^{n+1},y(t^{n+1}))]}
= y ( t n + 1 ) − y ( t n ) − h [ θ y ′ ( t n ) + ( 1 − θ ) y ′ ( t n + 1 ) ] {\displaystyle =y(t^{n+1})-y(t^{n})-h[\theta y'(t^{n})+(1-\theta )y'(t^{n+1})]}
= [ y ( t n ) + h y ′ ( t n ) + 1 2 h 2 y ″ ( t n ) + 1 6 h 3 y ′ ′ ′ ( t n ) ] − y ( t n ) − h { θ y ′ ( t n ) + ( 1 − θ ) [ y ′ ( t n ) + h y ″ ( t n ) + 1 2 h 2 y ′ ′ ′ ( t n ) ] } + O ( h 4 ) {\displaystyle =[y(t^{n})+hy'(t^{n})+{\frac {1}{2}}h^{2}y''(t^{n})+{\frac {1}{6}}h^{3}y^{\prime \prime \prime }(t^{n})]-y(t^{n})-h\{\theta y'(t^{n})+(1-\theta )[y'(t^{n})+hy''(t^{n})+{\frac {1}{2}}h^{2}y^{\prime \prime \prime }(t^{n})]\}+{\mathcal {O}}(h^{4})} = ( θ − 1 2 ) h 2 y ″ ( t n ) + ( 1 2 θ − 1 3 ) h 3 y ′ ′ ′ ( t n ) + O ( h 4 ) . {\displaystyle =\left(\theta -{\frac {1}{2}}\right)h^{2}y''(t^{n})+\left({\frac {1}{2}}\theta -{\frac {1}{3}}\right)h^{3}y^{\prime \prime \prime }(t^{n})+{\mathcal {O}}(h^{4}).}
從以下表達式中減去最後一個表示式：
y n + 1 − y n − h [ θ f ( t n , y n ) + ( 1 − θ ) f ( t n + 1 , y n + 1 ) ] = 0 , {\displaystyle y^{n+1}-y^{n}-h[\theta f(t^{n},y^{n})+(1-\theta )f(t^{n+1},y^{n+1})]=0,}
當 h {\displaystyle h}
e n + 1 , h = e n , h + θ h [ f ( t n , y ( t n ) + e n , h ) − f ( t n , y ( t n ) ) ] + ( 1 − θ ) h [ f ( t n + 1 , y ( t n + 1 ) + e n + 1 , h ) − f ( t n + 1 , y ( t n + 1 ) ) ] − 1 12 h 3 y ′ ′ ′ ( t n ) + O ( h 4 ) {\displaystyle e^{n+1,h}=e^{n,h}+\theta h[f(t^{n},y(t^{n})+e^{n,h})-f(t^{n},y(t^{n}))]+(1-\theta )h[f(t^{n+1},y(t^{n+1})+e^{n+1,h})-f(t^{n+1},y(t^{n+1}))]-{\frac {1}{12}}h^{3}y^{\prime \prime \prime }(t^{n})+{\mathcal {O}}(h^{4})} θ = 1 2 {\displaystyle \theta ={\frac {1}{2}}}
e n + 1 , h = e n , h + θ h [ f ( t n , y ( t n ) + e n , h ) − f ( t n , y ( t n ) ) ] + ( 1 − θ ) h [ f ( t n + 1 , y ( t n + 1 ) + e n + 1 , h ) − f ( t n + 1 , y ( t n + 1 ) ) ] + ( θ − 1 2 ) h 2 y ″ ( t n ) + O ( h 3 ) {\displaystyle e^{n+1,h}=e^{n,h}+\theta h[f(t^{n},y(t^{n})+e^{n,h})-f(t^{n},y(t^{n}))]+(1-\theta )h[f(t^{n+1},y(t^{n+1})+e^{n+1,h})-f(t^{n+1},y(t^{n+1}))]+(\theta -{\frac {1}{2}})h^{2}y''(t^{n})+{\mathcal {O}}(h^{3})} θ ≠ 1 2 {\displaystyle \theta \neq {\frac {1}{2}}}
其中 e i = y i − y ( t i ) {\displaystyle e^{i}=y^{i}-y(t^{i})} f {\displaystyle f} c {\displaystyle c} λ {\displaystyle \lambda }
‖ e n + 1 , h ‖ ≤ ‖ e n , h ‖ + θ h λ ‖ e n , h ‖ + ( 1 − θ ) h λ ‖ e n + 1 , h ‖ + c h 3 {\displaystyle \left\Vert e^{n+1,h}\right\Vert \leq \left\Vert e^{n,h}\right\Vert +\theta h\lambda \left\Vert e^{n,h}\right\Vert +(1-\theta )h\lambda \left\Vert e^{n+1,h}\right\Vert +ch^{3}} θ = 1 2 {\displaystyle \theta ={\frac {1}{2}}}
‖ e n + 1 , h ‖ ≤ ‖ e n , h ‖ + θ h λ ‖ e n , h ‖ + ( 1 − θ ) h λ ‖ e n + 1 , h ‖ + c h 2 {\displaystyle \left\Vert e^{n+1,h}\right\Vert \leq \left\Vert e^{n,h}\right\Vert +\theta h\lambda \left\Vert e^{n,h}\right\Vert +(1-\theta )h\lambda \left\Vert e^{n+1,h}\right\Vert +ch^{2}} θ ≠ 1 2 . {\displaystyle \theta \neq {\frac {1}{2}}.}
當 θ = 1 2 {\displaystyle \theta ={\frac {1}{2}}} [ 9] θ ≠ 1 2 {\displaystyle \theta \neq {\frac {1}{2}}} ‖ e n + 1 , h ‖ ≤ 1 + θ h λ 1 − ( 1 − θ ) h λ ‖ e n , h ‖ + c 1 − ( 1 − θ ) h λ h 2 . {\displaystyle \left\Vert e^{n+1,h}\right\Vert \leq {\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\left\Vert e^{n,h}\right\Vert +{\frac {c}{1-(1-\theta )h\lambda }}h^{2}.}
我們斷言 ‖ e n , h ‖ ≤ c λ [ ( 1 + θ h λ 1 − ( 1 − θ ) h λ ) n − 1 ] h {\displaystyle \left\Vert e^{n,h}\right\Vert \leq {\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{n}-1\right]h} n = 0 {\displaystyle n=0} ‖ e n , h ‖ = 0 {\displaystyle \left\Vert e^{n,h}\right\Vert =0} n = k {\displaystyle n=k} k ≥ 0 {\displaystyle k\geq 0} n = k + 1 {\displaystyle n=k+1}
‖ e k + 1 , h ‖ ≤ 1 + θ h λ 1 − ( 1 − θ ) h λ ‖ e k , h ‖ + c 1 − ( 1 − θ ) h λ h 2 , {\displaystyle \left\Vert e^{k+1,h}\right\Vert \leq {\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\left\Vert e^{k,h}\right\Vert +{\frac {c}{1-(1-\theta )h\lambda }}h^{2},}
然後代入
‖ e k n , h ‖ ≤ c λ [ ( 1 + θ h λ 1 − ( 1 − θ ) h λ ) k − 1 ] h . {\displaystyle \left\Vert e^{kn,h}\right\Vert \leq {\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{k}-1\right]h.}
我們有
‖ e k + 1 , h ‖ ≤ c λ [ ( 1 + θ h λ 1 − ( 1 − θ ) h λ ) k + 1 − 1 + θ h λ 1 − ( 1 − θ ) h λ ] h + c 1 − ( 1 − θ ) h λ h 2 {\displaystyle \left\Vert e^{k+1,h}\right\Vert \leq {\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{k+1}-{\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right]h+{\frac {c}{1-(1-\theta )h\lambda }}h^{2}}
= c λ [ ( 1 + θ h λ 1 − ( 1 − θ ) h λ ) k + 1 − 1 ] h . {\displaystyle ={\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{k+1}-1\right]h.}
所以我們的斷言對所有 n {\displaystyle n}
1 + θ h λ 1 − ( 1 − θ ) h λ = 1 + h λ 1 − ( 1 − θ ) h λ ≤ exp ( h λ 1 − ( 1 − θ ) h λ ) {\displaystyle {\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}{}=1+{\frac {h\lambda }{1-(1-\theta )h\lambda }}\leq \exp \left({\frac {h\lambda }{1-(1-\theta )h\lambda }}\right)}
根據指數函式的泰勒展開，我們有
‖ e n , h ‖ ≤ c λ [ ( 1 + θ h λ 1 − ( 1 − θ ) h λ ) n − 1 ] h {\displaystyle \left\Vert e^{n,h}\right\Vert {}\leq {\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{n}-1\right]h}
≤ c λ ( 1 + θ h λ 1 − ( 1 − θ ) h λ ) n h {\displaystyle \leq {\frac {c}{\lambda }}\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{n}h}
≤ c h λ exp ( n h λ 1 − ( 1 − θ ) h λ ) . {\displaystyle \leq {\frac {ch}{\lambda }}\exp \left({\frac {nh\lambda }{1-(1-\theta )h\lambda }}\right).}
根據我們的條件， n h ≤ t ∗ . {\displaystyle nh\leq t^{*}.} ‖ e n , h ‖ ≤ c h λ exp ( t ∗ λ 1 − ( 1 − θ ) h λ ) . {\displaystyle \left\Vert e^{n,h}\right\Vert \leq {\frac {ch}{\lambda }}\exp \left({\frac {t^{*}\lambda }{1-(1-\theta )h\lambda }}\right).} lim h → 0 ‖ e n , h ‖ = 0 {\displaystyle \lim _{h\rightarrow 0}\left\Vert e^{n,h}\right\Vert =0} 0 ≤ n h ≤ t ∗ {\displaystyle 0\leq nh\leq t^{*}} θ ≠ 1 2 . {\displaystyle \theta \neq {\frac {1}{2}}.}
需要注意的是，當 θ = 0 {\displaystyle \theta =0}
定理 後向尤拉法是一階收斂的。
備註 如果 f {\displaystyle f} f {\displaystyle f}
備註 如果不知道皮卡定理，可以透過時間離散化推匯出常微分方程解的存在性和唯一性。
現在我們考慮 y ′ = y 2 {\displaystyle y'=y^{2}} t ∈ [ 0 , 0.99 ] {\displaystyle t\in [0,0.99]} y ( t ) = 1 1 − t {\displaystyle y(t)={\frac {1}{1-t}}} 1 ≤ y ≤ 100 {\displaystyle 1\leq y\leq 100} f = y 2 {\displaystyle f=y^{2}} f {\displaystyle f} f ( x ) {\displaystyle f(x)} λ {\displaystyle \lambda } x , y ∈ R {\displaystyle x,y\in \mathbb {R} } ‖ f ( x ) − f ( y ) ‖ ‖ x − y ‖ ≤ λ {\displaystyle {\frac {\left\Vert f(x)-f(y)\right\Vert }{\left\Vert x-y\right\Vert }}\leq \lambda } x = 0 {\displaystyle x=0} ‖ y ‖ → ∞ , {\displaystyle \left\Vert y\right\Vert \rightarrow \infty ,} ‖ f ( x ) − f ( y ) ‖ ‖ x − y ‖ → ∞ . {\displaystyle {\frac {\left\Vert f(x)-f(y)\right\Vert }{\left\Vert x-y\right\Vert }}\rightarrow \infty .}
現在我們討論如何找到區域性 Lipschitz 常數 λ {\displaystyle \lambda } f {\displaystyle f} f {\displaystyle f} f {\displaystyle f}
‖ f ( y 1 ) − f ( y 2 ) ‖ = ‖ y 1 + y 2 ‖ ‖ y 1 − y 2 ‖ ≤ ( ‖ y 1 ‖ + ‖ y 2 ‖ ) ‖ y 1 − y 2 ‖ . {\displaystyle \left\Vert f(y^{1})-f(y^{2})\right\Vert =\left\Vert y^{1}+y^{2}\right\Vert \left\Vert y^{1}-y^{2}\right\Vert \leq \left(\left\Vert y^{1}\right\Vert +\left\Vert y^{2}\right\Vert \right)\left\Vert y^{1}-y^{2}\right\Vert .}
因此，找到 ‖ y ‖ {\displaystyle \left\Vert y\right\Vert } y ( t ) {\displaystyle y(t)} y ( t ) {\displaystyle y(t)} y ′ {\displaystyle y'} y ( t ) = 1 1 − t {\displaystyle y(t)={\frac {1}{1-t}}} 1 ≤ y ≤ 100 {\displaystyle 1\leq y\leq 100}
最後，我們得到了函式迭代和我們問題的後向尤拉方法的收斂性。因此，我們針對 y ′ = y 2 {\displaystyle y'=y^{2}} y ( 0 ) = 1 {\displaystyle y(0)=1} t ∈ [ 0 , 0.99 ] {\displaystyle t\in [0,0.99]}
推論 根據常微分方程解的存在性和唯一性定理以及本頁中的定理，我們得出最終目標，即後向尤拉方法與函式迭代生成的數值解在時間步長 h 0 {\displaystyle h0}
備註 這需要在進行函式迭代時仔細選擇初始迭代值。
備註 通常，ODE 的精確解是未知的，儘管可以推斷出區域性 Lipschitz 連續性。如果解變為無窮大，數值方法將無法收斂，或者如果近似解接近精確解，則會顯示非常大的值。在這些情況下解釋此類數值模擬時，需要謹慎。
以下兩個 Matlab 和 Python 程式演示了示例方程的後向尤拉法。第一個程式使用不動點迭代來求解非線性項，第二個程式使用牛頓法來求解非線性項。
一個演示不動點迭代的 Matlab 程式
程式碼下載
% A program to solve y'=y^2 using the backward Euler
% method and fixed point iteration
% This is not optimized and is very simple
clear all ; format compact ; format short ;
set ( 0 , 'defaultaxesfontsize' , 25 , 'defaultaxeslinewidth' , .7 , ...
'defaultlinelinewidth' , 6 , 'defaultpatchlinewidth' , 3.7 , ...
'defaultaxesfontweight' , 'bold' )
n = 10000 ; % Number of timesteps
Tmax = 0.99 ; % Maximum time
y0 = 1 ; % Initial value
tol = 0.1 ^10 ; % Tolerance for fixed point iterations
dt = Tmax / n ; % Time step
y = zeros ( 1 , n ); % vector for discrete solution
t = zeros ( 1 , n ); % vector for times of discrete solution
y ( 1 )= y0 ;
t ( 1 )= 0 ;
tic , % start timing
for i = 1 : n
yold = y ( i ); ynew = y ( i ); err = 1 ;
while err > tol
ynew = dt * yold ^2 + y ( i );
err = abs ( ynew - yold );
yold = ynew ;
end
y ( i + 1 )= ynew ;
t ( i + 1 )= t ( i ) + dt ;
end
toc , % stop timing
yexact = 1 ./ ( 1 - t ); max ( abs ( y - yexact )), % print the maximum error
figure ( 1 ); clf ; plot ( t , y , 'r+' , t , yexact , 'b-.' );
xlabel Time ; ylabel Solution ; legend ( 'Backward Euler' , 'Exact solution' );
title ( 'Numerical solution of dy/dt=y^2' );
一個演示不動點迭代的 Python 程式
程式碼下載
#!/usr/bin/env python
"""
A program to solve y'=y^2 using the backward Euler
method and fixed point iteration
This is not optimized and is very simple
"""
import time
import matplotlib.pyplot as plt
N = 1000 # Number of timesteps
tmax = 0.99 # Maximum time
y0 = 1
t0 = 0 # Initial value
tol = pow ( 0.1 , 10 ) # Tolerance for fixed point iterations
h = tmax / N # Time step
y = [ y0 ] # Variables holding the values of iterations
t = [ t0 ] # Times of discrete solutions
T0 = time . clock ()
for i in xrange ( N ):
yold = y [ i ]
ynew = y [ i ]
err = 1
while err > tol :
ynew = h * pow ( yold , 2 ) + y [ i ]
err = abs ( ynew - yold )
yold = ynew
y . append ( ynew )
t . append ( t [ i ] + h )
T = time . clock () - T0
yexact = [ 1.0 / ( 1.0 - x ) for x in t ]
print
print "Exact value: y( %d )= %f " % ( tmax , 1 / ( 1 - tmax ))
print "Value given by aproximation: y( %d )= %f " % ( tmax , y [ - 1 ])
maxerr = ( max ([ abs ( y [ i ] - yexact [ i ]) for i in xrange ( len ( y ))]))
print "Maximum error: %f " % maxerr
print "Elapsed time is %f " % ( T )
plt . figure ()
plt . plot ( t , y , 'r+' , t , yexact , 'b-.' )
plt . xlabel ( 'Time' )
plt . ylabel ( 'Solution' )
plt . legend (( 'Backward Euler' , 'Exact solution' ))
plt . title ( 'Numerical solution of dy/dt=y^2' )
plt . show ()
一個演示牛頓迭代的 Matlab 程式
程式碼下載
% A program to solve y'=y^2 using the backward Euler
% method and Newton iteration
% This is not optimized and is very simple
set ( 0 , 'defaultaxesfontsize' , 25 , 'defaultaxeslinewidth' , .7 , ...
'defaultlinelinewidth' , 6 , 'defaultpatchlinewidth' , 3.7 , ...
'defaultaxesfontweight' , 'bold' )
n = 100000 ; % Number of timesteps
Tmax = 0.99 ; % Maximum time
y0 = 1 ; % Initial value
tol = 0.1 ^10 ; % Tolerance for fixed point iterations
dt = Tmax / n ; % Time step
y = zeros ( 1 , n ); % vector for discrete solution
t = zeros ( 1 , n ); % vector for times of discrete solution
y ( 1 )= y0 ;
t ( 1 )= 0 ;
tic , % start timing
for i = 1 : n
yold = y ( i ); ynew = y ( i ); err = 1 ;
while err > tol
ynew = yold - ( yold - y ( i ) - dt * yold ^2 ) / ( 1 - 2 * dt * yold );
err = abs ( ynew - yold );
yold = ynew ;
end
y ( i + 1 )= ynew ;
t ( i + 1 )= t ( i ) + dt ;
end
toc , % stop timing
yexact = 1 ./ ( 1 - t ); max ( abs ( y - yexact )), % print maximum error
figure ( 1 ); clf ; plot ( t , y , 'r+' , t , yexact , 'b-.' );
xlabel Time ; ylabel Solution ;
legend ( 'Backward Euler' , 'Exact solution' );
title ( 'Numerical solution of dy/dt=y^2' );
一個演示牛頓迭代的 Python 程式。
程式碼下載
#!/usr/bin/env python
"""
A program to solve y'=y^2 using the backward Euler
method and Newton iteration
This is not optimized and is very simple
"""
import time
import matplotlib.pyplot as plt
N = 1000 # Number of timesteps
tmax = 0.99 # Maximum value
t0 = 0 # Initial t value
y0 = 1 # Initial value y(t0) = y0
tol = 0.1 ** 10 # Tolerance for fixed point iterations
h = ( tmax - t0 ) / N # Time step
y = [ y0 ] # List for discrete solutions
t = [ t0 ] # List with grid of discrete values of t
T0 = time . clock () #Start timing
for i in xrange ( N ):
yold = y [ i ]
ynew = y [ i ]
err = 1
while err > tol :
ynew = yold - ( yold - y [ i ] - h * ( yold ** 2 )) / ( 1 - 2 * h * yold )
err = abs ( ynew - yold )
yold = ynew
y . append ( ynew )
t . append ( t [ - 1 ] + h )
T = time . clock () - T0 # Stop timing
print "Exact value y( %f ) = %f " % ( t [ N ], 1 / ( 1 - t [ N ]))
print "Value given by approx y_n( %f ) = %f " % ( t [ N ], y [ N ])
print "The error = y-y_n = %f " % ( abs ( 1 / ( 1 - t [ N ]) - y [ N ]))
print "Time taken = %f " % ( T )
yexact = [ 1.0 / ( 1.0 - x ) for x in t ]
plt . figure ();
plt . plot ( t , y , 'r+' , t , yexact , 'b-.' );
plt . xlabel ( 'Time' )
plt . ylabel ( 'Solution' )
plt . legend (( 'Backward Euler' , 'Exact Solution' ))
plt . title ( 'Numerical solution of dy/dt=y^2' )
plt . show ()
執行 Matlab 中的不動點迭代程式，並檢查結果是否合理。現在，研究將時間步數從 0 到 0.99 的變化以及不動點迭代的容差對最大誤差的影響。特別是，嘗試 1,000 到 1,000,000（以 10 的冪表示）的時間步數範圍，以及 10 − 1 − 10 − 7 {\displaystyle 10^{-1}-10^{-7}} 10 − 1 {\displaystyle 10^{-1}}
編寫一個 Matlab 程式，使用 Crank-Nicolson 方法和不動點迭代來求解 y ′ = y 2 {\displaystyle y'=y^{2}} y ( 0 ) = 1 {\displaystyle y(0)=1}
證明微分方程 y ′ = | y | {\displaystyle y'={\sqrt {|y|}}} y ( 0 ) = 0 {\displaystyle y(0)=0}
找到該微分方程的至少兩個解析解。
使用前向尤拉法計算該微分方程的數值解。
使用後向尤拉法計算該微分方程的數值解。請務必嘗試不動點迭代的不同初始猜測，而不僅僅是前一步的值；您應該能夠計算出初始迭代的選擇對數值方法選擇的解的影響。對此發表評論。
使用隱式中點規則計算該微分方程的數值解。請務必嘗試不動點迭代的不同初始猜測，而不僅僅是前一步的值；您應該能夠計算出初始迭代的選擇對數值方法選擇的“解”的影響。對此發表評論。
使用牛頓迭代重複 (d) 和 (e)。
評論數值方法對求解沒有唯一解的微分方程的適用性。 修改一維 Allen-Cahn 方程的程式，使其使用 Crank-Nicolson 方法和非線性項的定點迭代。您需要在實空間中計算非線性項，以便您的最終方案為 u ^ n + 1 , k + 1 − u ^ n δ t = u ^ x x n + 1 , k + 1 + u ^ x x n 2 + 1 2 [ u n + 1 , k − ( u n + 1 , k ) 3 ] ^ + 1 2 [ u n − ( u n ) 3 ] ^ , {\displaystyle {\frac {{\hat {u}}^{n+1,k+1}-{\hat {u}}^{n}}{\delta t}}={\frac {{\hat {u}}_{xx}^{n+1,k+1}+{\hat {u}}_{xx}^{n}}{2}}+{\frac {1}{2}}{\widehat {\left[u^{n+1,k}-\left(u^{n+1,k}\right)^{3}\right]}}+{\frac {1}{2}}{\widehat {\left[u^{n}-\left(u^{n}\right)^{3}\right]}},} n {\displaystyle n} k {\displaystyle k}
修改二維 Allen-Cahn 方程的程式，使其使用 Crank-Nicolson 方法和非線性項的定點迭代。您需要在實空間中計算非線性項。
↑ Iserles (2009) ↑ Lax, Burstein 和 Lax (1979) ↑ Hughes 等人 (2008) ↑ Iserles (2009) ↑ Iserles (2009) ↑ Birkhoff 和 Rota (1989) ↑ Bradie (2006) ↑ Iserles (2009) ↑ Iserles (2009)
Birkhoff, G; Rota, G.C. (1989). 常微分方程 (第 4 版). 約翰·威立。
Bradie, B. (2006). 數值分析友好入門 . 皮爾森。
Hughes-Hallett, D. (2008). 微積分，單變數和多變數 (第 5 版). 約翰·威立。
Iserles, A. (2009). 微分方程數值分析入門 . 劍橋大學出版社。
Lax, A. (1976). 應用和計算微積分 . 第 1 卷. 施普林格。
![{\displaystyle t\in [t_{0},t^{*}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0aa5c94606c874bf4eaad4f5cb91bd9ee7471c23)
![{\displaystyle t\in [0,0.99]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1a02f3d61c8d8dde9b4457ba4023315143e77576)
![{\displaystyle [a,b]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9c4b788fc5c637e26ee98b45f89a5c08c85f7935)
![{\displaystyle x\in [a,b]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/026357b404ee584c475579fb2302a4e9881b8cce)
![{\displaystyle \left[x^{*}-\left|x^{*}-x_{0}\right|,x^{*}+\left|x^{*}-x_{0}\right|\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f9819f97dbda5be8675b2057d9095df7f16d0f26)
![{\displaystyle z_{n}\in [x_{n},x^{*}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6a7c51b6e7a0474705caf783053687259f20096f)
![{\displaystyle y^{n+1}=y^{n}+h[\theta f(t^{n},y^{n})+(1-\theta )f(t^{n+1},y^{n+1})]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3153cf734e644d4d57e298bc393bd1ade850fbef)
![{\displaystyle \theta \in [0,1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fead1e7dceab4be5ab2e91f5108144722daa8c36)
![{\displaystyle {}y(t^{n+1})-y(t^{n})-h[\theta f(t^{n},y(t^{n}))+(1-\theta )f(t^{n+1},y(t^{n+1}))]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/709823c215b17ad5bf0d74f446193ff2c2862f57)
![{\displaystyle =y(t^{n+1})-y(t^{n})-h[\theta y'(t^{n})+(1-\theta )y'(t^{n+1})]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d362c063380596d5cb80f6a2e0c40c241ff25eea)
![{\displaystyle =[y(t^{n})+hy'(t^{n})+{\frac {1}{2}}h^{2}y''(t^{n})+{\frac {1}{6}}h^{3}y^{\prime \prime \prime }(t^{n})]-y(t^{n})-h\{\theta y'(t^{n})+(1-\theta )[y'(t^{n})+hy''(t^{n})+{\frac {1}{2}}h^{2}y^{\prime \prime \prime }(t^{n})]\}+{\mathcal {O}}(h^{4})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e1bd05427d82f7bb1f082413ace5ca093b4197b1)
![{\displaystyle y^{n+1}-y^{n}-h[\theta f(t^{n},y^{n})+(1-\theta )f(t^{n+1},y^{n+1})]=0,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/32466900308f4c1a657f4b355cedd94e7cce94c4)
![{\displaystyle e^{n+1,h}=e^{n,h}+\theta h[f(t^{n},y(t^{n})+e^{n,h})-f(t^{n},y(t^{n}))]+(1-\theta )h[f(t^{n+1},y(t^{n+1})+e^{n+1,h})-f(t^{n+1},y(t^{n+1}))]-{\frac {1}{12}}h^{3}y^{\prime \prime \prime }(t^{n})+{\mathcal {O}}(h^{4})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fece5a3f9e49ebd0b59f8ccd9d68bee4bc777d43)
![{\displaystyle e^{n+1,h}=e^{n,h}+\theta h[f(t^{n},y(t^{n})+e^{n,h})-f(t^{n},y(t^{n}))]+(1-\theta )h[f(t^{n+1},y(t^{n+1})+e^{n+1,h})-f(t^{n+1},y(t^{n+1}))]+(\theta -{\frac {1}{2}})h^{2}y''(t^{n})+{\mathcal {O}}(h^{3})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/31243ccad0ac271d027644b28bbe0cc123ee2254)
![{\displaystyle \left\Vert e^{n,h}\right\Vert \leq {\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{n}-1\right]h}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4f6fc9636b8dc4f104c44b736fd5f32e5eea1e8c)
![{\displaystyle \left\Vert e^{kn,h}\right\Vert \leq {\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{k}-1\right]h.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ed341a28a88727035520be391291ef6b6dcd3a20)
![{\displaystyle \left\Vert e^{k+1,h}\right\Vert \leq {\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{k+1}-{\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right]h+{\frac {c}{1-(1-\theta )h\lambda }}h^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/01ba862390e928c6ee5e5db49099be3d42664991)
![{\displaystyle ={\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{k+1}-1\right]h.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8f8dfe84bc8d4643ca67f2809a4f9050da133dbb)
![{\displaystyle \left\Vert e^{n,h}\right\Vert {}\leq {\frac {c}{\lambda }}\left[\left({\frac {1+\theta h\lambda }{1-(1-\theta )h\lambda }}\right)^{n}-1\right]h}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b4db29f66d3da9e30bd38793847989c102438245)
![{\displaystyle {\frac {{\hat {u}}^{n+1,k+1}-{\hat {u}}^{n}}{\delta t}}={\frac {{\hat {u}}_{xx}^{n+1,k+1}+{\hat {u}}_{xx}^{n}}{2}}+{\frac {1}{2}}{\widehat {\left[u^{n+1,k}-\left(u^{n+1,k}\right)^{3}\right]}}+{\frac {1}{2}}{\widehat {\left[u^{n}-\left(u^{n}\right)^{3}\right]}},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e2929cf931d6af02269c314f54a8467d692fb13b)
