数值分析复习笔记-第九章-常微分方程初边值问题

Chapter9 常微分方程初边值问题

9.1 前言

1. 对微分方程，主要有两类问题：
   • 初值问题
     $$\{\begin{array}{r}{y}^{\prime }=f(x,y)\\ y(a)=y_0\end{array}$$
   • 边界值问题
     $$\{\begin{array}{rl}{y}^{\prime }& =f(x,y)\\ y(a)& =\alpha ,y(b)=\beta \end{array}$$
2. 其中，$x\in [a,b]$，$y\in R^m$
   $$\{\begin{array}{rl}m& =1,微分方程\\ m& >1,微分方程组\end{array}$$
3. 本章节考点：
   • 欧拉公式及其改进：欧拉公式、局部截断、p阶方法
   • 龙格-库塔公式：泰勒级数、二阶Runge-kutta
   • 单步法的收敛性和稳定性：会用即可

9.2 欧拉法及其改进

9.2.1 欧拉法

1. 很简单，就是一阶差分+近似：
   $$\begin{gathered} y^{\prime }(x_n)=f(x_n,y(x_n)) \\ \downarrow \\ h=x_{n+1}-x_n=\frac{b-a}{N} \\ \downarrow \\ {y}^{\prime }(x_n)\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})\approx y(x_n)+hf(x_n,y(x_n)) \\ \downarrow \\ y(x_{n+1})\approx y_{n+1} \\ \downarrow \\ {y}_{n+1}=y_n+hf(x_n,y_n) \end{gathered}$$
2. 几何意义：化曲为直->用多段折线来反映曲线
3. 实例代码：
   • 给定整数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 精度

1. 局部截断误差：假设前n次精确，第n+1次的误差=>$\widehat{y_{n+1}}$是用前n次真解算的
   $$y(x_{n+1})-\widehat{y_{n+1}}={\epsilon }_{n+1}$$
2. 整体截断误差：前n次累积误差，$y_n$是第n次近似解，误差会逐层累积
   $$y(x_n)-y_n=e_n$$
3. p-阶方法：局部截断误差为${\epsilon }_{n+1}=O(h^{p+1})$，即局部截断误差是$h^{p+1}$的等价无穷小，这个就称为p阶精度和p阶方法
4. 常用方法：

   欧拉方法：1阶方法
   $${\epsilon }_{n+1}=y(x_{n+1})-\widehat{y_{n+1}}$$
   $$y(x_{n+1})=y(x_n)+h{y}^{\prime }(x_{n+1})+\frac{h^2}{2}{y}^{\prime \prime }(\xi )$$
   $$\widehat{y_{n+1}}=y(x_n)+hf(x_n,y(x_n))=y(x_n)+h{y}^{\prime }(x_{n+1})$$
   $$|{\epsilon }_{n+1}|=|y(x_{n+1})-\widehat{y_{n+1}}|=\frac{h^2}{2}|y^{\prime \prime }(\xi )|\le \frac{M}{2}{h}^2=O(h^2)$$
   后退欧拉方法：2阶方法

9.2.3 欧拉法改进

1. 后退欧拉法：隐式公式->右端隐含$y_{n+1}$
   $$y_{n+1}=y_n+hf(x_{n+1},y_{n+1})$$
2. 梯形欧拉法：也是隐式公式
   $$y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}))$$
3. 改进的欧拉方程：显式公式，精度介于欧拉和后退欧拉公式，但是速度也介于欧拉和后退欧拉
   • 预估：$\widehat{y_{n+1}}=y_n+f(x_n,y_n)$
   • 校正：$y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},\widehat{y_{n+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})$在$x=x_n$处泰勒展开：
$$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 )$$

1. 舍去余项；
2. 对各项用已知条件表达；
   • p=1时，欧拉方法–>
     $$y\left(x_{n+1}\right)=y\left(x_n\right)+h{y}^{\prime }\left(x_n\right)$$
   • p=2时，对$f(x,y)$在$x$和$y$的偏导为$f_x,f_y$
     $$y^{\prime }(x_n)=f_n$$
     $$y^{\prime \prime }(x_n)=\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial x}=f_x+f_y{f}_n$$

9.3.2 基本思想：找平均斜率$K^{\ast }$

对公式：
$$\frac{y(x_{n+1})-y(x_n)}{h}=y^{\prime }(x_n+\theta h)=K^{\ast }$$
龙格库塔公式的基本思想就是找这个平均斜率$K^{\ast }$

1. 设法计算$f(x,y)$在某些节点的函数值，然后对这些函数值做线性组合，构造待定参数的近似计算公式
2. 然后将这些计算公式和真解泰勒展开比较，确定参数的取值，从而获得一定精度的计算公式

9.3.3 二段Runge-Kutta方法构造

1. p阶r段龙格-库塔公式：计算r次$f(x,y)$的值，且局部截断误差达到$O(h^{p+1})$
2. 二段龙格-库塔公式：
   $$\{\begin{array}{rl}{y}_{n+1}=& y_n+{\lambda }_1{k}_1+{\lambda }_2{k}_2\\ k_1=& hf\left(x_n,y_n\right)\\ k_2=& hf\left(x_n+{\alpha }_{2}h,y_n+{\beta }_{21}{k}_1\right)\end{array}$$
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_{n+1}=y_n+{\lambda }_1{k}_1+{\lambda }_2{k}_2=y_n+{\lambda }_{1}hf\left(x_n,y_n\right)+{\lambda }_{2}hf\left(x_n+{\alpha }_{2}h,y_n+{\beta }_{21}{k}_1\right)$$
   • 对$f\left(x_n+{\alpha }_{2}h,y_n+{\beta }_{21}{k}_1\right)$一阶泰勒展开：(对x和y分别求导展开)
     $$f\left(x_n+{\alpha }_{2}h,y_n+{\beta }_{21}{k}_1\right)=f_n+\alpha h{f}_x+{\beta }_{21}h{f}_y{f}_n+O(h^2)$$
   • 假设$y_n=y(x_n)$，化简近似公式：
     $$y_{n+1}=y(x_n)+({\lambda }_1+{\lambda }_2)f_{n}h+{\lambda }_2(\alpha f_x+{\beta }_{21}{f}_y{f}_n)h^2+O(h^3)$$
   • 对比：
     $$\begin{array}{rl}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_{n}h+{\lambda }_2(\alpha f_x+{\beta }_{21}{f}_y{f}_n)h^2+O(h^3)\end{array}$$
   • 得到：
     $$\{\begin{array}{c}{\lambda }_1+{\lambda }_2=1\\ {\alpha }_2={\beta }_{21}=\frac{1}{2{\lambda }_2}\end{array}$$

9.3.4 特点

1. 函数值次数增加，不能通过选择其他参数，来获得更高阶的龙格库塔公式
2. 可以通过改变${\lambda }_{1}and{\lambda }_2$取值来获得不同的段阶龙格库塔公式
3. 标准四阶四段龙格-库塔公式
   $$y_{n+1}=y_n+\frac{h}{6}\left(k_1+2k_2+2k_3+k_4\right)$$
   $$\begin{array}{rlrl}{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{array}$$

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 单步法

定义：

1. 单步法：已知$y_n$，推出$y_{n+1}$，此为单步法
2. 多步法：除了知道$y_n$，还要知道$y_{n-1}$、$y_{n-2}$等，才能推出$y_{n+1}$，此为多步法

9.4.2 收敛性

1. 定义：当步长$h\to 0$时，$y_n\to y(x_n)$
2. 定理：

对单步法：
$$y_{n+1}=y_n+h\phi \left(x_n,y_n,h\right)$$
$\phi$ 为增量函数，它依赖于f且仅是关于x、y、h的连续函数\

假设单步法是p阶函数，且增量函数满足李普希茨条件
$$|\phi (x,y,h)-\phi (x,\bar{y},h)|\le L_{\phi }|y-\bar{y}|$$
则有如下结论：
  (1) 该数值方法总体截断误差为$e_n=y\left(x_n\right)-y_n=O\left(h^p\right)$
  (2) $P\ge 1$，则该数值方法收敛

3. 改善数值方法收敛性
   • 加大N，即加密网格
   • 加阶级，提高数值方法的阶级数p

9.4.3 稳定性

1. 定义：用数值方法求解问题时，舍入误差会累积，当计算时的舍入误差被控制，即为稳定
   $${\delta }_m=|y_n-\widehat{y_n}|\le |{\delta }_n|$$
2. 影响因素：方法、步长、右端项$f(x,y)$
3. 绝对稳定域：在一个复平面内，数值方法稳定，在平面外不稳定
   • 对以下方程讨论稳定性：
     $$y^{\prime }=\lambda y(Re\lambda <0)$$
   • 欧拉方法绝对稳定域：$|1+h\lambda |\le 1$ -> 复平面以-1为圆心，1为半径的圆
   • 后退欧拉公式绝对稳定域：$|h\lambda -1|\ge 1$ -> 复平面以1为圆心，1为半径的外圆
4. A稳定性：某数值方法绝对稳定域包含$\lambda h$左半平面，则称为A稳定
   • 欧拉方法不具有A-稳定性
   • 后退欧拉方法具有A-稳定性
5. 由此可知，隐方法远远比显方法稳定