Chapter9 常微分方程初边值问题
9.1 前言
- 对微分方程,主要有两类问题:
- 初值问题
{ y ′ = f ( x , y ) y ( a ) = y 0 \left \{ \begin{aligned} y'=f(x,y) \\ y(a)=y_0 \end{aligned} \right. {y′=f(x,y)y(a)=y0 - 边界值问题
{ y ′ = f ( x , y ) y ( a ) = α , y ( b ) = β \left \{ \begin{aligned} y'&=f(x,y) \\ y(a)&=\alpha,y(b)=\beta \end{aligned} \right. {y′y(a)=f(x,y)=α,y(b)=β
- 初值问题
- 其中,
x
∈
[
a
,
b
]
x\in[a,b]
x∈[a,b],
y
∈
R
m
y\in R^m
y∈Rm
{ m = 1 , 微分方程 m > 1 , 微分方程组 \left \{ \begin{aligned} m&=1,微分方程 \\ m&>1,微分方程组 \end{aligned} \right. {mm=1,微分方程>1,微分方程组 - 本章节考点:
- 欧拉公式及其改进:欧拉公式、局部截断、p阶方法
- 龙格-库塔公式:泰勒级数、二阶Runge-kutta
- 单步法的收敛性和稳定性:会用即可
9.2 欧拉法及其改进
9.2.1 欧拉法
- 很简单,就是一阶差分+近似:
y ′ ( x n ) = f ( x n , y ( x n ) ) ↓ h = x n + 1 − x n = b − a N ↓ y ′ ( x n ) ≈ y ( x n + 1 ) − y ( x n ) x n + 1 − x n = y ( x n + 1 ) − y ( x n ) h ↓ y ( x n + 1 ) ≈ y ( x n ) + h f ( x n , y ( x n ) ) ↓ y ( x n + 1 ) ≈ y n + 1 ↓ y n + 1 = y n + h f ( x n , y n ) y'(x_n)=f(x_n,y(x_n))\\ \downarrow\\ h=x_{n+1}-x_n=\frac{b-a}{N}\\ \downarrow\\ y'(x_n)\boldsymbol{\approx}\frac{y(x_{n+1})-y(x_{n})}{x_{n+1}-x_n}=\frac{y(x_{n+1})-y(x_{n})}{h}\\ \downarrow\\ y(x_{n+1})\boldsymbol{\approx}y(x_{n})+hf(x_n,y(x_n))\\ \downarrow\\ y(x_{n+1})\boldsymbol{\approx}y_{n+1}\\ \downarrow\\ y_{n+1}=y_{n}+hf(x_n,y_{n}) y′(xn)=f(xn,y(xn))↓h=xn+1−xn=Nb−a↓y′(xn)≈xn+1−xny(xn+1)−y(xn)=hy(xn+1)−y(xn)↓y(xn+1)≈y(xn)+hf(xn,y(xn))↓y(xn+1)≈yn+1↓yn+1=yn+hf(xn,yn) - 几何意义:化曲为直->用多段折线来反映曲线
- 实例代码:
- 给定整数N:用来计算步长h=(b-a)/N
- 遍历x:计算y
function [x,y]=odeEuler(f,y0,a,b,N) h=(b-a)/N; x=a:h:b; y(1)=y0; for n=1:N y(n+1)=y(n)+h*feval(f,x(n),y(n)); end
9.2.2 精度
- 局部截断误差:假设前n次精确,第n+1次的误差=>
y
n
+
1
^
\widehat{y_{n+1}}
yn+1
是用前n次真解算的
y ( x n + 1 ) − y n + 1 ^ = ε n + 1 y(x_{n+1})-\widehat{y_{n+1}}=\varepsilon_{n+1} y(xn+1)−yn+1 =εn+1 - 整体截断误差:前n次累积误差,
y
n
y_{n}
yn是第n次近似解,误差会逐层累积
y ( x n ) − y n = e n y(x_{n})-y_{n}=e_{n} y(xn)−yn=en - p-阶方法:局部截断误差为 ε n + 1 = O ( h p + 1 ) \varepsilon_{n+1}=O(h^{p+1}) εn+1=O(hp+1),即局部截断误差是 h p + 1 h^{p+1} hp+1的等价无穷小,这个就称为p阶精度和p阶方法
- 常用方法:
欧拉方法:1阶方法
ε n + 1 = y ( x n + 1 ) − y n + 1 ^ \varepsilon_{n+1}=y(x_{n+1})-\widehat{y_{n+1}} εn+1=y(xn+1)−yn+1
y ( x n + 1 ) = y ( x n ) + h y ′ ( x n + 1 ) + h 2 2 y ′ ′ ( ξ ) y(x_{n+1})=y(x_{n})+hy'(x_{n+1})+\frac{h^2}{2}y''(\xi) y(xn+1)=y(xn)+hy′(xn+1)+2h2y′′(ξ)
y n + 1 ^ = y ( x n ) + h f ( x n , y ( x n ) ) = y ( x n ) + h y ′ ( x n + 1 ) \widehat{y_{n+1}}=y(x_{n})+hf(x_n,y(x_n))=y(x_{n})+hy'(x_{n+1}) yn+1 =y(xn)+hf(xn,y(xn))=y(xn)+hy′(xn+1)
∣ ε n + 1 ∣ = ∣ y ( x n + 1 ) − y n + 1 ^ ∣ = h 2 2 ∣ y ′ ′ ( ξ ) ∣ ≤ M 2 h 2 = O ( h 2 ) |\varepsilon_{n+1}|=|y(x_{n+1})-\widehat{y_{n+1}}|=\frac{h^2}{2}|y''(\xi)|\leq\frac{M}{2}h^2=O(h^2) ∣εn+1∣=∣y(xn+1)−yn+1 ∣=2h2∣y′′(ξ)∣≤2Mh2=O(h2)
后退欧拉方法:2阶方法
9.2.3 欧拉法改进
- 后退欧拉法:隐式公式->右端隐含
y
n
+
1
y_{n+1}
yn+1
y n + 1 = y n + h f ( x n + 1 , y n + 1 ) y_{n+1}=y_{n}+hf(x_{n+1},y_{n+1}) yn+1=yn+hf(xn+1,yn+1) - 梯形欧拉法:也是隐式公式
y n + 1 = y n + h 2 ( f ( x n , y n ) + f ( x n + 1 , y n + 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)) - 改进的欧拉方程:显式公式,精度介于欧拉和后退欧拉公式,但是速度也介于欧拉和后退欧拉
- 预估: y n + 1 ^ = y n + f ( x n , y n ) \widehat{y_{n+1}}=y_n+f(x_n,y_n) yn+1 =yn+f(xn,yn)
- 校正: y n + 1 = y n + h 2 ( f ( x n , y n ) + f ( x n + 1 , y n + 1 ^ ) ) y_{n+1}=y_{n}+\frac{h}{2}(f(x_{n},y_{n})+f(x_{n+1},\widehat{y_{n+1}})) yn+1=yn+2h(f(xn,yn)+f(xn+1,yn+1 ))
- 代码:预估步骤+校正步骤
function [x,y]=odegEuler(f,y0,a,b,N) h=(b-a)/N; x=a:h:b; y(1)=y0; for n=1:N yp=y(n)+h*feval(f,x(n),y(n)); yc=y(n)+h*feval(f,x(n),yp); y(n+1)=0.5*(yp+yc); end
9.3 Runge-Kutta法
9.3.1 Tarlor级数
y
(
x
n
+
1
)
y(x_{n+1})
y(xn+1)在
x
=
x
n
x=x_n
x=xn处泰勒展开:
y
(
x
n
+
1
)
=
y
(
x
n
)
+
h
y
′
(
x
n
)
+
h
2
2
y
′
′
(
x
n
)
+
⋯
+
h
p
p
!
y
(
p
)
(
x
n
)
+
h
p
+
1
(
p
+
1
)
!
y
(
p
+
1
)
(
ξ
)
y\left(x_{n+1}\right)=y\left(x_n\right)+h y^{\prime}\left(x_n\right)+\frac{h^2}{2} y^{\prime \prime}\left(x_n\right)+\cdots+\frac{h^p}{p !} y^{(p)}\left(x_n\right)+\frac{h^{p+1}}{(p+1) !} y^{(p+1)}(\xi)
y(xn+1)=y(xn)+hy′(xn)+2h2y′′(xn)+⋯+p!hpy(p)(xn)+(p+1)!hp+1y(p+1)(ξ)
- 舍去余项;
- 对各项用已知条件表达;
- p=1时,欧拉方法–>
y ( x n + 1 ) = y ( x n ) + h y ′ ( x n ) y\left(x_{n+1}\right)=y\left(x_n\right)+h y^{\prime}\left(x_n\right) y(xn+1)=y(xn)+hy′(xn) - p=2时,对
f
(
x
,
y
)
f(x,y)
f(x,y)在
x
x
x和
y
y
y的偏导为
f
x
,
f
y
f_x,f_y
fx,fy
y ′ ( x n ) = f n y'(x_n)=f_n y′(xn)=fn
y ′ ′ ( x n ) = ∂ f ∂ x + ∂ f ∂ y ∂ y ∂ x = f x + f y f n y''(x_n)=\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial x}=f_x+f_yf_n y′′(xn)=∂x∂f+∂y∂f∂x∂y=fx+fyfn
- p=1时,欧拉方法–>
9.3.2 基本思想:找平均斜率 K ∗ K^* K∗
对公式:
y
(
x
n
+
1
)
−
y
(
x
n
)
h
=
y
′
(
x
n
+
θ
h
)
=
K
∗
\frac{y(x_{n+1})-y(x_n)}{h}=y'(x_n+\theta h)=K^*
hy(xn+1)−y(xn)=y′(xn+θh)=K∗
龙格库塔公式的基本思想就是找这个平均斜率
K
∗
K^*
K∗
- 设法计算 f ( x , y ) f(x,y) f(x,y)在某些节点的函数值,然后对这些函数值做线性组合,构造待定参数的近似计算公式
- 然后将这些计算公式和真解泰勒展开比较,确定参数的取值,从而获得一定精度的计算公式
9.3.3 二段Runge-Kutta方法构造
- p阶r段龙格-库塔公式:计算r次 f ( x , y ) f(x,y) f(x,y)的值,且局部截断误差达到 O ( h p + 1 ) O(h^{p+1}) O(hp+1)
- 二段龙格-库塔公式:
{ y n + 1 = y n + λ 1 k 1 + λ 2 k 2 k 1 = h f ( x n , y n ) k 2 = h f ( x n + α 2 h , y n + β 21 k 1 ) \left\{\begin{aligned} y_{n+1}=&y_n+\lambda_1 k_1+\lambda_2 k_2 \\ k_1 =&h f\left(x_n, y_n\right) \\ k_2 =&h f\left(x_n+\alpha_2 h, y_n+\beta_{21} k_1\right) \end{aligned}\right. ⎩ ⎨ ⎧yn+1=k1=k2=yn+λ1k1+λ2k2hf(xn,yn)hf(xn+α2h,yn+β21k1) - 推导如下:
- 真值泰勒展开:
y ( x n + 1 ) = y ( x n ) + h y ′ ( x n ) + h 2 2 y ′ ′ ( x n ) + O ( h 3 ) = y ( x n ) + h f n + h 2 2 ( f x + f y f ) + O ( h 3 ) y\left(x_{n+1}\right)=y\left(x_n\right)+h y^{\prime}\left(x_n\right)+\frac{h^2}{2} y^{\prime \prime}\left(x_n\right)+O(h^{3})=y(x_n)+h f_n+\frac{h^2}{2}\left(f_x+f_y f\right)+O(h^{3}) y(xn+1)=y(xn)+hy′(xn)+2h2y′′(xn)+O(h3)=y(xn)+hfn+2h2(fx+fyf)+O(h3) - 对龙格-库塔公式合并:
y n + 1 = y n + λ 1 k 1 + λ 2 k 2 = y n + λ 1 h f ( x n , y n ) + λ 2 h f ( x n + α 2 h , y n + β 21 k 1 ) y_{n+1}=y_n+\lambda_1 k_1+\lambda_2 k_2=y_n+\lambda_1 h f\left(x_n, y_n\right)+\lambda_2 h f\left(x_n+\alpha_2 h, y_n+\beta_{21} k_1\right) yn+1=yn+λ1k1+λ2k2=yn+λ1hf(xn,yn)+λ2hf(xn+α2h,yn+β21k1) - 对
f
(
x
n
+
α
2
h
,
y
n
+
β
21
k
1
)
f\left(x_n+\alpha_2 h, y_n+\beta_{21} k_1\right)
f(xn+α2h,yn+β21k1)一阶泰勒展开:(对x和y分别求导展开)
f ( x n + α 2 h , y n + β 21 k 1 ) = f n + α h f x + β 21 h f y f n + O ( h 2 ) f\left(x_n+\alpha_2 h, y_n+\beta_{21} k_1\right)=f_n+\alpha hf_x+\beta_{21} hf_yf_n+O(h^{2}) f(xn+α2h,yn+β21k1)=fn+αhfx+β21hfyfn+O(h2) - 假设
y
n
=
y
(
x
n
)
y_n=y(x_n)
yn=y(xn),化简近似公式:
y n + 1 = y ( x n ) + ( λ 1 + λ 2 ) f n h + λ 2 ( α f x + β 21 f y f n ) h 2 + O ( h 3 ) y_{n+1}=y(x_n)+(\lambda_1+\lambda_2) f_nh+ \lambda_2(\alpha f_x+\beta_{21} f_yf_n)h^2+O(h^{3}) yn+1=y(xn)+(λ1+λ2)fnh+λ2(αfx+β21fyfn)h2+O(h3) - 对比:
y ( x n + 1 ) = y n + h f n + h 2 2 ( f x + f y f ) + O ( h 3 ) y n + 1 = y ( x n ) + ( λ 1 + λ 2 ) f n h + λ 2 ( α f x + β 21 f y f n ) h 2 + O ( h 3 ) \begin{aligned} y\left(x_{n+1}\right)&=y_n+h f_n+\frac{h^2}{2}\left(f_x+f_y f\right)+O(h^{3})\\ y_{n+1}&=y(x_n)+(\lambda_1+\lambda_2) f_nh+ \lambda_2(\alpha f_x+\beta_{21} f_yf_n)h^2+O(h^{3}) \end{aligned} y(xn+1)yn+1=yn+hfn+2h2(fx+fyf)+O(h3)=y(xn)+(λ1+λ2)fnh+λ2(αfx+β21fyfn)h2+O(h3) - 得到:
{ λ 1 + λ 2 = 1 α 2 = β 21 = 1 2 λ 2 \left\{\begin{array}{c} \lambda_1+\lambda_2=1 \\ \alpha_2=\beta_{21}=\frac{1}{2 \lambda_2} \end{array}\right. {λ1+λ2=1α2=β21=2λ21
- 真值泰勒展开:
9.3.4 特点
- 函数值次数增加,不能通过选择其他参数,来获得更高阶的龙格库塔公式
- 可以通过改变 λ 1 a n d λ 2 \lambda_1 and \lambda_2 λ1andλ2取值来获得不同的段阶龙格库塔公式
- 标准四阶四段龙格-库塔公式
y n + 1 = y n + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) y_{n+1}=y_n+\frac{h}{6}\left(k_1+2 k_2+2 k_3+k_4\right) yn+1=yn+6h(k1+2k2+2k3+k4)
k 1 = f ( y n , t n ) k 2 = f ( y n + h k 1 2 , t n + h 2 ) k 3 = f ( y n + h k 2 2 , t n + h 2 ) k 4 = f ( y n + h k 3 , t n + h ) \begin{aligned} k_1 & =f\left(y_n, t_n\right) & k_2 & =f\left(y_n+h \frac{k_1}{2}, t_n+\frac{h}{2}\right) \\ k_3 & =f\left(y_n+h \frac{k_2}{2}, t_n+\frac{h}{2}\right) & k_4 & =f\left(y_n+h k_3, t_n+h\right) \end{aligned} k1k3=f(yn,tn)=f(yn+h2k2,tn+2h)k2k4=f(yn+h2k1,tn+2h)=f(yn+hk3,tn+h)
9.3.5 四阶-四段龙格-库塔公式代码
function [x,y]=rk4(f,y0,a,b,N)
h=(b-a)/N;
x=a:h:b;
y(1)=y0;
for n=1:N
k1=h*feval(f,x(n),y(n));
k2=h*feval(f,x(n)+0.5*h,y(n)+0.5*k1);
k3=h*feval(f,x(n)+0.5*h,y(n)+0.5*k2);
k4=h*feval(f,x(n)+h,y(n)+k3);
y(n+1)=y(n)+(k1+k2+k3+k4)/6;
end
9.4 单步法的收敛性和稳定性
9.4.1 单步法
定义:
- 单步法:已知 y n y_n yn,推出 y n + 1 y_{n+1} yn+1,此为单步法
- 多步法:除了知道 y n y_{n} yn,还要知道 y n − 1 y_{n-1} yn−1、 y n − 2 y_{n-2} yn−2等,才能推出 y n + 1 y_{n+1} yn+1,此为多步法
9.4.2 收敛性
- 定义:当步长 h → 0 h→0 h→0时, y n → y ( x n ) y_{n}\rightarrow y(x_n) yn→y(xn)
- 定理:
对单步法:
y n + 1 = y n + h φ ( x n , y n , h ) y_{n+1}=y_n+h \varphi\left(x_n, y_n, h\right) yn+1=yn+hφ(xn,yn,h)
φ \varphi φ 为增量函数,它依赖于f且仅是关于x、y、h的连续函数\假设单步法是p阶函数,且增量函数满足李普希茨条件
∣ φ ( 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ˉ∣
则有如下结论:
(1) 该数值方法总体截断误差为 e n = y ( x n ) − y n = O ( h p ) e_n=y\left(x_n\right)-y_n=O\left(h^p\right) en=y(xn)−yn=O(hp)
(2) P ≥ 1 P\geq1 P≥1,则该数值方法收敛
- 改善数值方法收敛性
- 加大N,即加密网格
- 加阶级,提高数值方法的阶级数p
9.4.3 稳定性
- 定义:用数值方法求解问题时,舍入误差会累积,当计算时的舍入误差被控制,即为稳定
δ m = ∣ y n − y n ^ ∣ ≤ ∣ δ n ∣ \delta_m=|y_n-\widehat{y_n}|\leq |\delta_n| δm=∣yn−yn ∣≤∣δn∣ - 影响因素:方法、步长、右端项 f ( x , y ) f(x,y) f(x,y)
- 绝对稳定域:在一个复平面内,数值方法稳定,在平面外不稳定
- 对以下方程讨论稳定性:
y ′ = λ y ( R e λ < 0 ) y'=\lambda y\space \space(Re \lambda<0) y′=λy (Reλ<0) - 欧拉方法绝对稳定域: ∣ 1 + h λ ∣ ≤ 1 |1+h\lambda|\leq1 ∣1+hλ∣≤1 -> 复平面以-1为圆心,1为半径的圆
- 后退欧拉公式绝对稳定域: ∣ h λ − 1 ∣ ≥ 1 |h\lambda -1| \geq 1 ∣hλ−1∣≥1 -> 复平面以1为圆心,1为半径的外圆
- 对以下方程讨论稳定性:
- A稳定性:某数值方法绝对稳定域包含
λ
h
\lambda h
λh左半平面,则称为A稳定
- 欧拉方法不具有A-稳定性
- 后退欧拉方法具有A-稳定性
- 由此可知,隐方法远远比显方法稳定