隐式显式方法避免了直接求解非线性问题。这对某些问题可能是有利的，但在其他问题中可能导致嚴重的步长限制。此外，由此产生的数值方案有时可能具有不良的定性特性。出于这个原因，我们需要描述允许我们求解完全隐式数值方案中产生的非线性方程的方法。
我们考虑一个常微分方程 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)
