\qquad常微分方程是描述连续变化的数学语言,微分方程的求解就是确定给定方程的可微方程y(t)y(t)y(t)。考虑一阶常微分方程的初值问题
y′=f(x,y),x∈[x0,b],(1)y(x0)=y0(2)
y^{'}=f(x,y),\quad x\in [x_{0},b],\quad(1)\\
y(x_{0})=y_{0}\quad(2)
y′=f(x,y),x∈[x0,b],(1)y(x0)=y0(2)\qquad如果存在实数L>0L>0L>0,使得
∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣,∀y1,y2∈R,
|f(x,y_{1})-f(x,y_{2})|\leq L|y_{1}-y_{2}|,\forall y_{1},y_{2}\in\mathbb{R},
∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣,∀y1,y2∈R,则称fff关于yyy满足利普西茨(LipschiteLipschiteLipschite)条件,LLL称为fff的利普西茨常数(简称Lisp.Lisp.Lisp.常数)
\qquad定理1 设fff在区域D={(x,y)∣a≤x≤b,y∈R}D=\{(x,y)|a\leq x\leq b,y\in\mathbb{R}\}D={(x,y)∣a≤x≤b,y∈R}上连续,关于yyy满足利普西茨条件,则对任意x0∈[a,b],y0∈Rx_{0}\in [a,b],y_{0}\in\mathbb{R}x0∈[a,b],y0∈R,常微分方程初值问题满足(1)(2)式当x∈[a,b]x\in [a,b]x∈[a,b]时存在唯一的连续可微解y(x)y(x)y(x)。
\qquad所谓数值解法,就是寻求解y(x)y(x)y(x)在一系列离散节点
x1<x2<⋯<xn<⋯
x_{1}<x_{2}<\cdots<x_{n}<\cdots
x1<x2<⋯<xn<⋯\quad上的近似值y1,y2,⋯ ,yn,⋯y_{1},y_{2},\cdots,y_{n},\cdotsy1,y2,⋯,yn,⋯相邻两节点的间距hn=xn+1−xnh_{n}=x_{n+1}-x_{n}hn=xn+1−xn称为步长。求解过程首先对常微分方程离散化,建立求数值解的递推公式。当计算yn+1y_{n+1}yn+1时只用到前一点yny_{n}yn,称为单步法。用到前面kkk个点的值yn,yn−1,⋯y_{n},y_{n-1},\cdotsyn,yn−1,⋯时称为kkk步法。
\qquad欧拉法 设做出折线(由切线相连(x0,x1,⋯xnx_{0},x_{1},\cdots x_{n}x0,x1,⋯xn处作切线))的顶点PnP_{n}Pn,过Pn(xn,yn)P_{n}(x_{n},y_{n})Pn(xn,yn)依方向场的方向推进到Pn+1(xn+1,yn+1)P_{n+1}(x_{n+1},y_{n+1})Pn+1(xn+1,yn+1),显然两个顶点Pn,Pn+1P_{n},P_{n+1}Pn,Pn+1的坐标有关系
yn+1−ynxn+1−xn=f(xn,yn)
\frac{y_{n+1}-y_{n}}{x_{n+1}-x_{n}}=f(x_{n},y_{n})
xn+1−xnyn+1−yn=f(xn,yn)\qquad对微分方程式(1)中的导数用均差近似,即
yn+1−ynh≈y′(xn)=f(xn,y(xn))
\frac{y_{n+1}-y_{n}}{h}\approx y^{'}(x_{n})=f(x_{n},y(x_{n}))
hyn+1−yn≈y′(xn)=f(xn,y(xn))\qquad依均差近似逐次计算。为分析计算公式的精度,使用泰勒展开将yn+1y_{n+1}yn+1在xnx_{n}xn处展开,则有
y(xn+1)=y(xn+h)=y(xn)+y′(xn)h+h22y′′(ξn),ξn∈(xn,xn+1)
y(x_{n+1})=y(x_{n}+h)=y(x_{n})+y^{'}(x_{n})h+\frac{h^{2}}{2}y^{''}(\xi_{n}),\quad \xi_{n}\in(x_{n},x_{n+1})
y(xn+1)=y(xn+h)=y(xn)+y′(xn)h+2h2y′′(ξn),ξn∈(xn,xn+1)\qquad在yn=yxny_{n}=y_{x_{n}}yn=yxn的前提下,f(xn,yn)=f(xn,y(xn))=y′(xn)f(x_{n},y_{n})=f(x_{n},y(x_{n}))=y^{'}(x_{n})f(xn,yn)=f(xn,y(xn))=y′(xn),于是欧拉法的误差
y(xn+1)−yn+1=h22y′′(ξn)≈h22y′′(xn)
y(x_{n+1})-y_{n+1}=\frac{h^{2}}{2}y^{''}(\xi_{n})\approx\frac{h^{2}}{2}y^{''}(x_{n})
y(xn+1)−yn+1=2h2y′′(ξn)≈2h2y′′(xn)称为此方法的局部截断误差。显式单步法的局部截断误差一般性的表示方法为:
Tn+1=y(xn+1)−y(xn)−hφ(xn,y(xn),h)=O(hp+1)
T_{n+1}=y(x_{n+1})-y(x_{n})-h\varphi(x_{n},y(x_{n}),h)=O(h^{p+1})
Tn+1=y(xn+1)−y(xn)−hφ(xn,y(xn),h)=O(hp+1)则称该方法具有ppp阶精度。
\qquad如果对微分方程从xnx_{n}xn到xn+1x_{n+1}xn+1积分,得
y(xn+1)=y(xn)+∫xnxn+1f(t,y(t))dt
y(x_{n+1})=y(x_{n})+\int_{x_{n}}^{x_{n+1}}f(t,y(t))dt
y(xn+1)=y(xn)+∫xnxn+1f(t,y(t))dt右端积分使用右矩形公式hf(xn+1,y(xn+1))hf(x_{n+1},y(x_{n+1}))hf(xn+1,y(xn+1))近似得到另一个公式
yn+1=yn+hf(xn+1,yn+1)
y_{n+1}=y_{n}+hf(x_{n+1},y_{n+1})
yn+1=yn+hf(xn+1,yn+1)通常使用迭代法求解(实现xn→xn+1x_{n}\to x_{n+1}xn→xn+1的运算),设用欧拉公式
yn+1(0)=yn+hf(xn,yn)
y_{n+1}^{(0)}=y_{n}+hf(x_{n},y_{n})
yn+1(0)=yn+hf(xn,yn)给出迭代初值yn+1(0)y_{n+1}^{(0)}yn+1(0)带入上式得
yn+1(1)=yn+hf(xn+1,yn+1(0))
y_{n+1}^{(1)}=y_{n}+hf(x_{n+1},y^{(0)}_{n+1})
yn+1(1)=yn+hf(xn+1,yn+1(0))反复迭代得到
yn+1(k+1)=yn+hf(xn+1,yn+1(k)),k=0,1,⋯
y_{n+1}^{(k+1)}=y_{n}+hf(x_{n+1},y^{(k)}_{n+1}),k=0,1,\cdots
yn+1(k+1)=yn+hf(xn+1,yn+1(k)),k=0,1,⋯由利普西茨条件得(证明迭代法收敛\color{#F00}{证明迭代法收敛}证明迭代法收敛)
∣yn+1(k+1)−yn+1∣=h∣f(xn+1,yn+1(k))−f(xn+1,yn+1)∣≤hL∣yn+1(k)−yn+1∣
|y^{(k+1)}_{n+1}-y_{n+1}|=h|f(x_{n+1},y^{(k)}_{n+1})-f(x_{n+1},y^{}_{n+1})|\leq hL|y^{(k)}_{n+1}-y_{n+1}|
∣yn+1(k+1)−yn+1∣=h∣f(xn+1,yn+1(k))−f(xn+1,yn+1)∣≤hL∣yn+1(k)−yn+1∣由此,只要hL<1hL<1hL<1使用迭代法就能收敛到解yn+1y_{n+1}yn+1。
\qquad梯形方法(增加精度) 公式为yn+1=yn+h2[f(xn,yn)+f(xn+1,yn+1)]y_{n+1} = y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},y_{n+1})]yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1)],迭代公式为:
yn+1(0)=yn+hf(xn,yn)yn+1(k+1)=yn+h2(f(xn,yn)+f(xn+1,yn+1(k)),k=0,1,2,⋯
y_{n+1}^{(0)}=y_{n}+hf(x_{n},y_{n})\\
y_{n+1}^{(k+1)}=y_{n}+\frac{h}{2}(f(x_{n},y_{n})+f(x_{n+1},y_{n+1}^{(k)}),\quad k=0,1,2,\cdots
yn+1(0)=yn+hf(xn,yn)yn+1(k+1)=yn+2h(f(xn,yn)+f(xn+1,yn+1(k)),k=0,1,2,⋯\qquad改进欧拉公式(增加了欧拉公式的精度,减小梯形算法迭代带来的计算次数),建立预测-校正系统。
预测 yˉn+1=yn+hf(xn,yn)\bar{y}_{n+1}=y_{n}+hf(x_{n},y_{n})yˉn+1=yn+hf(xn,yn)
校正 yn+1=yn+h2[f(xn,yn)+f(xn+1,yˉn+1)]y_{n+1}=y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},\bar{y}_{n+1})]yn+1=yn+2h[f(xn,yn)+f(xn+1,yˉn+1)]
\qquad龙格-库塔法,由公式
y(xn+1)=y(xn)+∫xnxn+1f(t,y(t))dt
y(x_{n+1})=y(x_{n})+\int_{x_{n}}^{x_{n+1}}f(t,y(t))dt
y(xn+1)=y(xn)+∫xnxn+1f(t,y(t))dt要使近似精度提高就是提高右式的近似精度,近似公式为:
∫xnxn+1f(t,y(t))dt≈h∑i=1rcif(xn+λih,y(xn+λih))
\int_{x_{n}}^{x_{n+1}}f(t,y(t))dt\approx h\sum_{i=1}^{r}c_{i}f(x_{n}+\lambda _{i}h,y(x_{n}+\lambda _{i}h))
∫xnxn+1f(t,y(t))dt≈hi=1∑rcif(xn+λih,y(xn+λih))分割的段越多,rrr越大,上式右端(相当于增量函数φ(xn,y(xn)h)\varphi(x_{n},y(x_{n})h)φ(xn,y(xn)h))的近似精度越高。由公式
yn+1=yn+hφ(xn,yn,h)
y_{n+1} = y_{n} + h\varphi(x_{n},y_{n},h)
yn+1=yn+hφ(xn,yn,h)其中,
φ(xn,yn,h)=∑i=1rciKiK1=f(xn,yn)Ki=f(xn+λih,yn+∑j=1i−1μijKj)i=2,3,⋯ ,r
\varphi(x_{n},y_{n},h)=\sum_{i=1}^{r}c_{i}K_{i}\\
K_{1}=f(x_{n},y_{n})\\
K_{i}=f(x_{n}+\lambda_{i}h,y_{n}+\sum_{j=1}^{i-1}\mu_{ij}K_{j})\quad i=2,3,\cdots,r
φ(xn,yn,h)=i=1∑rciKiK1=f(xn,yn)Ki=f(xn+λih,yn+j=1∑i−1μijKj)i=2,3,⋯,r 通过给定局部截断误差精度求解未知系数。对于步长的选择,由给定的精度ε\varepsilonε和偏差Δ=∣yn+1(h2)−yn+1(h)∣\Delta=|y_{n+1}^{(\frac{h}{2})}-y_{n+1}^{(h)}|Δ=∣yn+1(2h)−yn+1(h)∣(折半前后两次的偏差\color{#F00}{折半前后两次的偏差}折半前后两次的偏差)决定。当精度大于偏差时,减小步长(对折),反之则反。
\qquad单步法的收敛性 假设单步法具有ppp阶精度,且增量函数φ(x,y,h)\varphi(x,y,h)φ(x,y,h)关于yyy满足利普西茨条件
∣φ(x,y,h)−φ(x,yˉ,h)∣≤Lφ∣y−yˉ∣
|\varphi(x,y,h)-\varphi(x,\bar{y},h)|\leq L_{\varphi}|y-\bar{y}|
∣φ(x,y,h)−φ(x,yˉ,h)∣≤Lφ∣y−yˉ∣又设初值y0y_{0}y0是准确的,即y0=y(x0)y_{0}=y(x_{0})y0=y(x0),其整体的截断误差
y(xn)−yn=O(hp)
y(x_{n})-y_{n}=O(h^{p})
y(xn)−yn=O(hp)稳定性 由模型方程y′=λyy^{'}=\lambda yy′=λy的欧拉公式为:
yn+1=(1+hλ)yn
y_{n+1}=(1+h\lambda)y_{n}
yn+1=(1+hλ)yn设在节点值yny_{n}yn上有一个扰动值εn\varepsilon_{n}εn,它的传播使节点yn+1y_{n+1}yn+1产生εn+1\varepsilon_{n+1}εn+1的扰动值,扰动值满足
εn+1=(1+hλ)εn
\varepsilon_{n+1}=(1+h\lambda)\varepsilon_{n}
εn+1=(1+hλ)εn为使差分方程的解非增(解稳定),选择hhh充分小,使得∣1+hλ∣≤1|1+h\lambda|\leq1∣1+hλ∣≤1,即−2<λ<0-2<\lambda<0−2<λ<0。