7.1 - 7.4 差分公式

一阶常微分方程初值问题

一阶常微分方程初值问题的一般形式为：
\[ \left\{\begin{array}{l}\displaystyle\frac{dy}{dx}=f(x,y),\quad a\le x\le b\\y(a)=\alpha\end{array}\right. \] 其中 \(f(x,y)\) 为已知函数， \(\alpha\) 为给定的值。

在许多数学模型中，上述方程通常以 \(x\) 描述时间，而定解条件 \(y(a)=\alpha\) 则给出了函数 \(y(x)\) 在初始时刻的取值。因此称为初值问题。

问题： + 上述方程何时存在唯一解 + 如何计算 \(y(x)\)

Lipschitz条件： 若函数 \(f(x,y)\) 在区域 \(\{a\le x\le b,\ m<y<M\}\) 上连续，满足 \[ \forall y,\bar{y},\ |f(x,y)-f(x,\bar{y})|\le L|y-\bar{y}| \] 其中 \(L>0\) 为Lipschitz常数（此处Lipschitz常数可以 \(\ge1\) ），则初值问题在初始时刻 \(a\) 的某个邻域上存在唯一解。 （不满足Lipschitz条件时，不一定存在唯一解。）

构造一阶常微分方程初值问题数值解法

假设初值问题的解 \(y=y(x)\) 唯一存在且足够光滑。对求解区域 \([a,b]\) 做等距剖分 \(a=x_0<x_1<x_2<\dots<x_n<\dots<x_N=b\) 。称 \(h=(b-a)/N\) 为剖分步长， \(x_n=a+nh,\ n=0,1,\dots,N\) 为剖分节点。数值解法即求精确解 \(y(x)\) 在剖分节点 \(x_n\) 上的值 \(y(x_n)\) 的近似值 \(y_n\) 。

差分公式：在区间 \([x_n,x_{n+1}]\) 上对微分方程两端同时积分有： \[ \displaystyle y(x_{n+1})-y(x_n)=\int_{x_n}^{x_{n+1}}f(x,y(x))dx \] 对该式右边积分部分应用不同的数值积分公式（参考前一章）做逼近，就得到相应不同的差分公式。

  • Euler公式：对右边积分应用左矩形公式 \(\displaystyle\int_a^bf(x)dx\approx(b-a)f(a)\) ，得到Euler差分公式： \[ \left\{\begin{array}{l}y_{n+1}=y_n+hf(x_n,y_n)\\y_0=\alpha\\n=0,1,\dots,N-1\end{array}\right. \]
  • 梯形公式：对右边积分应用梯形公式 \(\displaystyle\int_a^bf(x)dx\approx\frac{(b-a)}{2}[f(a)+f(b)]\) ，得到梯形差分公式： \[ \left\{\begin{array}{l}\displaystyle y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_{n+1})]\\y_0=\alpha\\n=0,1,\dots,N-1\end{array}\right. \]
  • Euler中点公式：将积分范围扩大为 \([x_{n-1},x_{n+1}]\) ，有 \[ \displaystyle y(x_{n+1})-y(x_{n-1})=\int_{x_{n-1}}^{x_{n+1}}f(x,y(x))dx \] 对右边积分应用中矩形公式 \(\displaystyle\int_a^bf(x)dx\approx(b-a)[f(\frac{a+b}{2})]\) ，得到Euler中点公式（或称双步Euler公式）： \[ \left\{\begin{array}{l}\displaystyle y_{n+1}=y_{n-1}+2hf(x_n,y_n)\\y_0=\alpha\\n=0,1,\dots,N-1\end{array}\right. \] 该公式属于多步方法，需要更多初值信息。

其中，Euler公式和Euler中点公式为显式方法，梯形公式为隐式方法。

改进Euler方法

梯形公式计算精度好，但属于隐式公式，不便计算。

根据非线性方程迭代法的思想，进行如下近似计算： \[ \left\{\begin{aligned} \displaystyle&y_{n+1}^{[0]}=y_n+hf(x_n)\\ \displaystyle&y_{n+1}^{[k+1]}=y_n+\frac{h}{2}[f(x_n,y_n),f(x_{n+1},y_{n+1}^{[k]})]\\ &&k=0,1,\dots\\ &y_0=\alpha,\ n=0,1,\dots,N-1 \end{aligned}\right. \] 即首先应用Euler公式提供 \(y_{n+1}\) 的初始值（预估），然后采用梯形公式框架进行关于 \(y_{n+1}\) 值的迭代计算（校正）。

校正过程中迭代计算是否一定收敛？

考察迭代格式，迭代函数为 \[ \displaystyle\varphi(x)=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y)] \]

压缩映射原理：设函数 \(f\) 定义域和值域均为 \([a,b]\) ，并存在一个常数 \(L\) ，满足 \(0<L<1\) ，使得对 \(\forall x,y\in[a,b]\) ，都有 \(|f(x)-f(y)|\le L|x-y|\) ，则称 \(f\) 是 \([a,b]\) 上的一个压缩映射，称常数 \(L\) 为Lipschitz常数（压缩常数）。

假设 \(\displaystyle\frac{\partial f}{\partial y}\) 存在，则当 \(|\varphi'(y)|\le L<1\) 也即 \(\displaystyle\frac{h}{2}\left|\frac{\partial f}{\partial y}\right|\le L<1\) 时，迭代必然收敛。理论上只需令剖分步长 \(h\) 足够小即可满足条件。 计算中，当 \(\left|y_{n+1}^{[k+1]}-y_{n+1}^{[k]}\right|<\varepsilon\) ，取 \(y_{n+1}=y_{n+1}^{[k+1]}\)

若仅迭代一步则有： \[ \left\{\begin{aligned} \displaystyle&\bar{y}_{n+1}=y_n+hf(x_n)\\ \displaystyle&y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n),f(x_{n+1},\bar{y}_{n+1})]\\ &y_0=\alpha,\ n=0,1,\dots,N-1 \end{aligned}\right. \] 称之为改进Euler方法，属于单步显式方法。也可写作： \[ \left\{\begin{aligned} \displaystyle&y_{n+1}=y_n+h(K_1+K_2)\\ &K_1=f(x_n,y_n)\\ &K_2=f(x_n+h,y_n+hK_1)\\ &y_0=\alpha,\ n=0,1,\dots,N-1 \end{aligned}\right. \]

误差分析

在节点 \(x_{n+1}\) 处的误差 \(y(x_{n+1})-y_{n+1}\) ，不仅与 \(y_{n+1}\) 这一步计算，而且与前 \(n\) 步均有关。为简化误差分析，着重研究一步计算时产生的截断误差，假设 \(y_n=y(x_n)\) ，称 \(y(x_{n+1})-y_{n+1}\) 为局部截断误差。

若单步差分公式的局部截断误差为 \(O(h^{p+1})\) （同阶无穷小），则称该公式为 \(p\) 阶方法。 \(p\) 为非负整数，阶数越高精度越好。

由一元Taylor公式： \[ \displaystyle y(x_{n+1})=y(x_n+h)=y(x_n)+y'(x_n)h+\frac{y''(x_n)}{2!}h^2+\frac{y'''(x_n)}{3!}h^3+O(h^4) \] 由二元Taylor公式： \[ \displaystyle\begin{aligned}f(x_n+h,y_n+k)=&f(x_n,y_n)+\frac{\partial f(x_n,y_n)}{\partial x}h+\frac{\partial f(x_n,y_n)}{\partial y}k\\ &+\frac{1}{2!}\left[\frac{\partial^2f(x_n,y_n)}{\partial x^2}h^2+2\frac{\partial^2 f(x_n,y_n)}{\partial x \partial y}hk+\frac{\partial^2f(x_n,y_n)}{\partial y^2} k^2\right]\\ &+\dots+\frac{1}{k!}\left[h\frac{\partial}{\partial x}+k\frac{\partial}{\partial y}\right]^kf(x_n,y_n)+\dots \end{aligned} \]

常见差分公式的局部截断误差 + Euler公式： \[ \displaystyle y(x_{n+1})-y_{n+1}=\frac{y''(x_n)}{2!}h^2+O(h^3)=O(h^2) \] 因此Euler公式为 \(1\) 阶方法。 + 改进Euler公式： \[ y(x_{n+1})-y_{n+1}=O(h^3) \] 因此Euler公式为 \(2\) 阶方法。