miniyuan

常微分方程初值问题的数值解法

题 9.1

令 f(x,y)=y−2xyf(x, y) = y - \dfrac{2x}{y}，x0=0x_0 = 0，y0=1y_0 = 1。

  1. Euler 法：

    yn+1=yn+hf(xn,yn)y_{n+1} = y_n + h f(x_n, y_n)

    从而：

    y1=y0+hf(x0,y0)=1+0.1×1=1.1y_1 = y_0 + hf(x_0, y_0) = 1 + 0.1 \times 1 = 1.1

    解析解为：

    y(x1)=1.2≈1.095445y(x_1) = \sqrt{1.2} \approx 1.095445

    绝对误差：

    ∣e1∣=∣y(x1)−y1∣≈0.004555|e_1| = |y(x_1) - y_1| \approx 0.004555

  2. 改进 Euler 法

    yn+1=yn+h2[f(xn,yn)+f(xn+1,yn+hf(xn,yn))]y_{n+1} = y_n + \frac{h}{2}[f(x_n, y_n) + f(x_{n+1}, y_n + hf(x_n, y_n))]

    其中：

    y0+hf(x0,y0)=1.1y_0 + hf(x_0, y_0) = 1.1

    从而：

    y1=y0+h2[f(x0,y0)+f(x1,y0+hf(x0,y0))]=1+0.05×(1+f(0.1,1.1))=1+0.05×(1+101110)=24112200≈1.095909\begin{aligned} y_1 &= y_0 + \frac{h}{2}[f(x_0, y_0) + f(x_1, y_0 + hf(x_0, y_0))] \\ &= 1 + 0.05 \times (1 + f(0.1, 1.1)) \\ &= 1 + 0.05 \times (1 + \frac{101}{110}) \\ &= \frac{2411}{2200} \approx 1.095909 \end{aligned}

    解析解为：

    y(x1)=1.2≈1.095445y(x_1) = \sqrt{1.2} \approx 1.095445

    绝对误差：

    ∣e2∣=∣y(x1)−y1∣≈0.000464|e_2| = |y(x_1) - y_1| \approx 0.000464

    可见改进 Euler 法的精度比 Euler 法高一个数量级。

  3. 四阶 Runge-Kutta 法：

    k1=hf(xn,yn)k2=hf(xn+h2, yn+12k1)k3=hf(xn+h2, yn+12k2)k4=hf(xn+h, yn+k3)yn+1=yn+16(k1+2k2+2k3+k4)\begin{aligned} k_1 &= hf(x_n, y_n) \\ k_2 &= hf\left(x_n+\frac{h}{2},\ y_n+\frac{1}{2}k_1\right) \\ k_3 &= hf\left(x_n+\frac{h}{2},\ y_n+\frac{1}{2}k_2\right) \\ k_4 &= hf\left(x_n+h,\ y_n+k_3\right) \\ y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned}

    故：

    k1=0.2×1=0.2k2=0.2×f(0.1,1.1)=0.2×(1.1−0.21.1)≈0.183636k3=0.2×f(0.1,1.091818)≈0.181727k4=0.2×f(0.2,1.181727)≈0.168648\begin{aligned} k_1 &= 0.2 \times 1 = 0.2 \\ k_2 &= 0.2 \times f(0.1, 1.1) = 0.2 \times (1.1 - \frac{0.2}{1.1}) \approx 0.183636 \\ k_3 &= 0.2 \times f(0.1, 1.091818) \approx 0.181727 \\ k_4 &= 0.2 \times f(0.2, 1.181727) \approx 0.168648 \end{aligned}

    从而：

    y1=y0+16(k1+2k2+2k3+k4)≈1.183229y_1 = y_0 + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \approx 1.183229

    解析解为：

    y(x1)=1.4≈1.183216y(x_1) = \sqrt{1.4} \approx 1.183216

    绝对误差：

    ∣e3∣=∣y(x1)−y1∣≈0.000013|e_3| = |y(x_1) - y_1| \approx 0.000013

题 9.2

局部截断误差为：

ln+1=y(xn+1)−y~n+1l_{n+1}=y(x_{n+1})-\tilde y_{n+1}

其中 y~n+1\tilde y_{n+1} 满足：

y~n+1=y(xn)+h[c f(xn,y(xn))+(1−c) f(xn+1,y~n+1)]\tilde y_{n+1}=y(x_n)+h\Bigl[c\,f\bigl(x_n,y(x_n)\bigr)+(1-c)\,f\bigl(x_{n+1},\tilde y_{n+1}\bigr)\Bigr]

对 y(xn+1)y(x_{n+1}) 在 xnx_n 处展开：

y(xn+1)=y(xn)+hy′(xn)+h22y′′(xn)+h36y′′′(xn)+O(h4)y(x_{n+1})=y(x_n)+hy'(x_n)+\frac{h^2}{2}y''(x_n)+\frac{h^3}{6}y'''(x_n)+O(h^4)

对 y(xn)y(x_n) 在 xn+1x_{n+1} 处展开：

y(xn)=y(xn+1)−hy′(xn+1)+h22y′′(xn+1)−h36y′′′(xn+1)+O(h4)y(x_n)=y(x_{n+1})-hy'(x_{n+1})+\frac{h^2}{2}y''(x_{n+1})-\frac{h^3}{6}y'''(x_{n+1})+O(h^4)

也即：

y(xn+1)=y(xn)+hy′(xn+1)−h22y′′(xn+1)+h36y′′′(xn+1)+O(h4)y(x_{n+1})=y(x_n)+hy'(x_{n+1})-\frac{h^2}{2}y''(x_{n+1})+\frac{h^3}{6}y'''(x_{n+1})+O(h^4)

将两式分别乘以 cc 与 (1−c)(1-c) 后相加，得到 y(xn+1)y(x_{n+1}) 的表示：

y(xn+1)=y(xn)+h[cy′(xn)+(1−c)y′(xn+1)]+h22[cy′′(xn)−(1−c)y′′(xn+1)]+h36[cy′′′(xn)+(1−c)y′′′(xn+1)]+O(h4)\begin{aligned} y(x_{n+1})=y(x_n) &+h\Bigl[cy'(x_n)+(1-c)y'(x_{n+1})\Bigr] +\frac{h^2}{2}\Bigl[cy''(x_n)-(1-c)y''(x_{n+1})\Bigr] \\ &+\frac{h^3}{6}\Bigl[cy'''(x_n)+(1-c)y'''(x_{n+1})\Bigr] +O(h^4) \end{aligned}

令：

y∗=y(xn)+h[c f(xn,y(xn))+(1−c) f(xn+1,y(xn+1))]=y(xn)+h[c y′(xn)+(1−c) y′(xn+1)]\begin{aligned} y^* &=y(x_n)+h\Bigl[c\,f\bigl(x_n,y(x_n)\bigr)+(1-c)\,f\bigl(x_{n+1},y(x_{n+1})\bigr)\Bigr] \\ &=y(x_n)+h\Bigl[c\,y'(x_n)+(1-c)\,y'(x_{n+1})\Bigr] \end{aligned}

由 y~n+1\tilde y_{n+1} 的定义与 y∗y^* 的定义相减：

y~n+1−y∗=h(1−c)[f(xn+1,y~n+1)−f(xn+1,y(xn+1))]=h(1−c)fy(xn+1,η)(y~n+1−y(xn+1))=−h(1−c)fy(xn+1,η)ln+1\begin{aligned} \tilde y_{n+1}-y^* &=h(1-c)\Bigl[f\bigl(x_{n+1},\tilde y_{n+1}\bigr)-f\bigl(x_{n+1},y(x_{n+1})\bigr)\Bigr] \\ &=h(1-c)f_y\bigl(x_{n+1},\eta\bigr)\bigl(\tilde y_{n+1}-y(x_{n+1})\bigr) \\ &=-h(1-c)f_y\bigl(x_{n+1},\eta\bigr)l_{n+1} \end{aligned}

从而：

ln+1=y(xn+1)−y∗+y∗−y~n+1=y(xn+1)−y∗+h(1−c)fy(xn+1,η)ln+1\begin{aligned} l_{n+1} &= y(x_{n+1})-y^*+y^*-\tilde{y}_{n+1} \\ &= y(x_{n+1})-y^*+h(1-c)f_y\bigl(x_{n+1},\eta\bigr)l_{n+1} \end{aligned}

也即：

ln+1[1−h(1−c)fy]=y(xn+1)−y∗l_{n+1}\Bigl[1-h(1-c)f_y\Bigr]=y(x_{n+1})-y^*

从而：

ln+1=y(xn+1)−y∗1−h(1−c)fy=[y(xn+1)−y∗][1+h(1−c)fy+O(h2)]=h22[cy′′(xn)−(1−c)y′′(xn+1)]+O(h3)=h2(c−12)y′′(xn+1)+O(h3)\begin{aligned} l_{n+1} &=\frac{y(x_{n+1})-y^*}{1-h(1-c)f_y} \\ &=\bigl[y(x_{n+1})-y^*\bigr]\Bigl[1+h(1-c)f_y+O(h^2)\Bigr] \\ &=\frac{h^2}{2}\Bigl[cy''(x_n)-(1-c)y''(x_{n+1})\Bigr]+O(h^3) \\ &=h^2\Bigl(c-\frac{1}{2}\Bigr)y''(x_{n+1})+O(h^3) \end{aligned}

  • 当 c≠12c\neq\dfrac12 时：

    ln+1=O(h2)l_{n+1}=O(h^2)

  • 当 c=12c=\dfrac12 时：

    ln+1=O(h3)l_{n+1}=O(h^3)

由上可知，对任意参数 cc，均有 ln+1=O(h2)l_{n+1}=O(h^2)，因此对任意 cc，该数值格式都是相容的。

且当 c=12c=\dfrac12 时，局部截断误差阶最高，为 2 阶。此时即为梯形公式。