[计算机数值分析]四阶龙格-库塔经典格式解常微分方程的初值问题

Spring-_-Bear 的 CSDN 博客导航

---

式（1）是四阶龙格 - 库塔格式的一种常用经典格式：

$$\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4)\\ K_1=f(x_n,y_n)\\ K_2=f(x_{n+\frac{1}{2}},y_n+\frac{h}{2}{K}_1)\\ K_3=f(x_{n+\frac{1}{2}},y_n+\frac{h}{2}{K}_2)\\ K_4=f(x_{n+1},y_n+h{K}_3)\end{array}\text{\hspace{1em}(1)}$$

经典的龙格 - 库塔格式（1）每一步需要 4 次计算函数值 $f$，可以直接验证它具有四阶精度，不过其论证极其繁琐，此处从略。下图描述了四阶经典的龙格 - 库塔方法。

例：设取步长 h = 0.2，从 x = 0 到 x = 1 用四阶经典格式（1）解决以下常微分方程的初值问题。

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

解：这里四阶经典格式（1）中 K₁，K₂，K₃，K₄ 的具体形式是

$$K_1=y_n-\frac{2x_n}{y_n}$$

$$K_2=y_n+\frac{h}{2}{K}_1-\frac{2x_n+h}{y_n+\frac{h}{2}{K}_1}$$

$$K_3=y_n+\frac{h}{2}{K}_2-\frac{2x_n+h}{y_n+\frac{h}{2}{K}_2}$$

$$K_4=y_n+h{K}_3-\frac{2(x_n+h)}{y_n+h{K}_3}$$

下表中记录了计算结果，其中 $y(x_n)$ 仍表示准确解。

xn	yn	y(xn)
0.2	1.1832	1.1832
0.4	1.3417	1.3416
0.6	1.4833	1.4832
0.8	1.6125	1.6125
1.0	1.7321	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 K1 = f(x0, y0);
        double K2 = f(x0 + h / 2, y0 + h / 2 * K1);
        double K3 = f(x0 + h / 2, y0 + h / 2 * K2);
        double K4 = f(x1, y0 + h * K3);
        // 求新值
        double y1 = y0 + h / 6 * (K1 + 2 * K2 + 2 * K3 + K4);
        // 输出本次步进后的离散点数据
        cout << i << "\t" << setw(10) << x1 << "\t" << setw(10) << y1 << "\t" << setw(10) << f1(x1) << endl;

        x0 = x1;
        y0 = y1;
    }

    return 0;
}