[计算机数值分析]改进欧拉格式解常微分方程的初值问题

Spring-_-Bear 的 CSDN 博客导航

---

欧拉方法

$$y_{n+1}=y_n+hf(x_n,y_n)，n=0,1,2,...$$

是一种显式算法，计算量小，但精度很低。梯形方法

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

虽然提高了精度，但它是一种隐式算法，需要借助于迭代过程求解，计算量大。

为此，我们将综合应用这两种方法，先用欧拉格式求得一个初步的近似值，记为 ${\overline{y}}_{n+1}$，称之为预报值；预报值 ${\overline{y}}_{n+1}$ 的精度不高，我们用它替代梯形公式

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

右端的 ${\overline{y}}_{n+1}$ 再直接计算，得到校正值 $y_{n+1}$。这样建立的预报 - 校正系统

$$\{\begin{array}{l}预报：{\overline{y}}_{n+1}=y_n+hf(x_n,y_n)\\ 校正：y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},{\overline{y}}_{n+1})]\end{array}$$

称作改进的欧拉格式。这是一种一步显式格式，它可表为如下嵌套格式

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

或表为下列平均化形式

$$\{\begin{array}{l}{y}_p=y_n+hf(x_n,y_n)\\ y_c=y_n+hf(x_{n+1},y_p)\\ y_{n+1}=\frac{1}{2}(y_p+y_c)\end{array}$$

图 3-2 描述了改进的欧拉方法，其中 h 为步长，N 为步数，x₀，y₀ 为 “老值”，即前一步的近似解，x₁，y₁ 为 “新值”，即该步计算结果。

例：用改进欧拉方法求解以下常微分方程的初值问题：

$$\{\begin{array}{l}{y}^{\prime }=y-\frac{2x}{y}(0<x<1)\\ y(0)=1\end{array}$$

解：求解初值问题的改进欧拉格式具有以下形式：

$$\{\begin{array}{l}{y}_p=y_n+h(y_n-\frac{2x_n}{y_n})\\ y_c=y_n+h(y_p-\frac{2x_{n+1}}{y_p})\\ y_{n+1}=\frac{1}{2}(y_p+y_c)\end{array}$$

取 h = 0.1，计算结果如下表所示，同欧拉方法的计算结果比较，改进欧拉方法明显地改善了精度。

xn	yn	y(xn)
0.1	1.0959	1.0954
0.2	1.1841	1.1832
0.3	1.2662	1.2649
0.4	1.3434	1.3416
0.5	1.4164	1.4142
0.6	1.4860	1.4832
0.7	1.5525	1.5492
0.8	1.6165	1.6125
0.9	1.6782	1.6733
1.0	1.7379	1.7321

运行示例：

程序源码：

#include <iostream>
#include <iomanip>
#include <math.h>

using namespace std;

/**
 * 欧拉格式中的 x，y 的函数关系式，即 f(xn,yn)
 */
double f(double x, double y)
{
    return y - 2 * x / y;
}

/**
 * 实际函数解析式
 */
double f1(double x)
{
    return sqrt(1 + 2 * x);
}

int main(void)
{
    double x0, y0;
    cout << "请输入初值:";
    cin >> x0 >> y0;

    double h;
    cout << "请输入步长:";
    cin >> h;

    int N;
    cout << "请输入步数:";
    cin >> N;

    cout << "\t" << setw(10) << "xn"
        << "\t" << setw(10) << "yn"
        << "\t" << setw(10) << "y(xn)" << endl;

    for (int i = 1; i <= N; i++)
    {
        // 离散点
        double x1 = x0 + h;
        // 根据欧拉格式求出的值，预报
        double yp = y0 + h * f(x0, y0);
        // 校正
        double yc = y0 + h * f(x1, yp);
        double y1 = (yp + yc) / 2;
        // 输出本次步进后的离散点数据
        cout << i << "\t" << setw(10) << x1 << "\t" << setw(10) << y1 << "\t" << setw(10) << f1(x1) << endl;

        x0 = x1;
        y0 = y1;
    }

    return 0;
}