GSL：求解具有时间相关系数的 ODE答案

【问题标题】：GSL: Solving ODEs with time-dependent coefficientsGSL：求解具有时间相关系数的 ODE
【发布时间】：2014-01-20 19:23:52
【问题描述】：

我有一个 ODE 类型：

x'(t) = a(t)x+g(t)

我要解决的问题。唯一的GSL ODE example 不是很有帮助，因为唯一的系数 (\mu) 与时间无关。

这个问题一直是answered on the GSL mailing list，但是答案很不清楚——g(t) 被忽略，也没有解释如何将a(t) 合并到func 中（是否应该在*params 中传递？）。

有没有我可以看到使用 GSL 解决此类 ODE 的示例？

更新： 正如下面所指出的，这已在 GSL 邮件列表中得到答复。这是一个完整的示例程序：

#include <stdio.h>
#include <math.h>
#include "gsl/gsl_errno.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_odeiv2.h"
int func(double t, const double y[], double f[], void *params) {
    f[0] = -t* y[0];
    return GSL_SUCCESS;
}
int jac(double t, const double y[], double *dfdy, double dfdt[], void
*params) {
    gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy, 1, 1);
    gsl_matrix * m = &dfdy_mat.matrix;
    gsl_matrix_set(m, 0, 0, -t);
    dfdt[0] = -1;
    return GSL_SUCCESS;
}
int main(void) {
    double mu = 0;
    gsl_odeiv2_system sys = { func, jac, 1, &mu };
    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new(&sys,
            gsl_odeiv2_step_rk1imp, 1e-7, 1e-7, 0.0);
    int i;
    double t = 0.0, t1 = 2.0;
    const double x0 = 1.0;
    double y[1] = {x0};
    const int N = 100;
    printf("time\t \tapprox solution \t exact solution\n");
    for (i = 0; i <= N; i++) {
        double ti = i * (t1 / N);
        int status = gsl_odeiv2_driver_apply(d, &t, ti, y);
        if (status != GSL_SUCCESS) {
            printf("error, return value=%d\n", status);
            break;
        }
        printf("%.5e\t%.5e\t\t%.5e\n", t, y[0], x0*exp(-0.5*t*t));
    }
    gsl_odeiv2_driver_free(d);
    printf("\n...and done!\n");
    return 0;
}

【问题讨论】：

标签： c++ gsl ode

【解决方案1】：

如果您不限于 GSL 和/或 C，您可以使用 http://odeint.com - 用于 ODE 的现代 C++ 库。 Odeint 是 boost 的一部分，因此它可能已经安装在您的系统上，或者可以轻松安装在大多数 Linux 发行版的包管理器中。

您可以简单地定义您的系数和 ODE，并使用例如 RK4 方法：

double coef_a( double t ) { /* return a(t) */ };
double coef_g( double t ) { /* return b(t) */ };

typedef std::array< double , 1 > state_type;
double ode( state_type const &x , state_type &dxdt , double )
{
    dxdt[0] = coef_a( t ) * x[0] + coef_g( t );
}

state_type x;
double t_state , t_end , dt;
// initialize x
integrate_const( runge_kutta< state_type >() , ode , x , t_start , dt , t_end );

【讨论】：

谢谢！对于不断发展的 ODE，GSL 与 ODEint 的优缺点是什么？
odeint 优点：现代 C++，高性能，适用于 GPU，非常灵活，odeint 缺点：难以阅读内部结构（如果您不习惯 C++），gsl 优点：更多隐式方法，gsl 缺点: 老式 C 风格的 API。当然，我的观点偏向于 odeint :)

【解决方案2】：

这里已经回答了这个问题http://lists.gnu.org/archive/html/help-gsl/2012-06/msg00036.html

【讨论】：

您应该将该问题的相关部分粘贴到您的答案中。
@Magnilex 是对的，我已将程序发布到我的问题中，但如果它在这里会更好。