基于GSL库的龙格库塔法解常微分方程实例

以下是一个基于GSL库的龙格库塔法解常微分方程的实例：

#include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

int func(double t, const double y[], double f[], void *params) {
    f[0] = y[1];
    f[1] = -y[0];
    return GSL_SUCCESS;
}

int main() {
    const double y0[] = { 1.0, 0.0 };
    const double t0 = 0.0;
    const double h = 0.1;
    const double tf = 10.0;

    gsl_odeiv2_system sys = { func, NULL, 2, NULL };

    gsl_odeiv2_driver *driver = gsl_odeiv2_driver_alloc_y_new(
        &sys, gsl_odeiv2_step_rk4, h, 1e-6, 1e-6);

    double t = t0;
    double y[2] = { y0[0], y0[1] };

    std::cout << "t\ty[0]\ty[1]" << std::endl;

    while (t < tf) {
        int status = gsl_odeiv2_driver_apply(driver, &t, t + h, y);

        if (status != GSL_SUCCESS) {
            std::cerr << "error: driver returned " << status << std::endl;
            break;
        }

        std::cout << t << "\t" << y[0] << "\t" << y[1] << std::endl;
    }

    gsl_odeiv2_driver_free(driver);
    return 0;
}

这个程序使用龙格库塔法（RK4）来解决常微分方程y'' + y = 0。函数func定义了ODE系统，其中f是输出函数，params是传递给函数的任意参数。在主函数中，我们定义了初始条件y0和时间步长h，然后使用gsl_odeiv2_driver来执行求解。在每个时间步骤中，我们将当前时间t和当前状态y传递给gsl_odeiv2_driver_apply，以计算下一个状态。