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

数值分析中，龙格－库塔法（Runge-Kutta methods）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。

Spring-_-Bear

5974人浏览 · 2021-12-14 10:34:03

Spring-_-Bear · 2021-12-14 10:34:03 发布

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

$\begin{cases} 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}+hK_{3}) \end{cases} \tag1$

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

在这里插入图片描述

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

$\begin{cases} y'=y-\frac{2x}{y} (0 < x < 1)\\ y(0)=1 \end{cases}$

解：这里四阶经典格式（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}+hK_{3}-\frac{2(x_{n}+h)}{y_{n}+hK_{3}}$

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

x_n	y_n	y(x_n)
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;
}