c - 如何查阅具有许多参数和谐波函数的 GSL ODE

Question

我正在使用 GSL 研究非线性微分方程。问题是我对 C 语言还很陌生。我刚刚将 GNU 网站上的示例改编成我现在感兴趣的方程。

这是等式:

d2x/dt2 + r*dx/dy + cos(x) + v*cos(2*x+0.4) E1*sin(wt) + E2*sin(2*w*t+a) = 0

我陷入困境的是我不知道如何在代码中插入多个参数。此外，我不知道如何在这段代码中使用余弦或正弦函数。

我试图通过Google一路搜索来解决这个问题。我找不到任何对我有帮助的东西。

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

int func (double t, const double x[], double y[], void *params)
{
    double r = *(double *)params;
    double v = *(double *)params;
    double w = *(double *)params;
    double E1 = *(double *)params;
    double E2  = *(double *)params;
    double a  = *(double *)params;
    y[0] = x[1];
    y[1] = -r*x[1] - cos(x[0]) - v*cos(2*x[0]+0.4) - E1*sin(w*t) - E2*sin(2*w*t+a);
    return GSL_SUCCESS;
}

int jac (double t, const double x[], double *dydx, double dydt[], void *params)
{
    double r = *(double *)params;
    double v = *(double *)params;
    double w = *(double *)params;
    double E1 = *(double *)params;
    double E2  = *(double *)params;
    double a  = *(double *)params;
    gsl_matrix_view dydx_mat = gsl_matrix_view_array (dydx, 2, 2);
    gsl_matrix * m = &dydx_mat.matrix;
    gsl_matrix_set (m, 0, 0, 0.0);
    gsl_matrix_set (m, 0, 1, 1.0);
    gsl_matrix_set (m, 1, 0, sin(x[0]) + 2*v*sin(2*x[0]+0.4));
    gsl_matrix_set (m, 1, 1, -r);
    dydt[0] = 0.0;
    dydt[1] = 0.0;
    return GSL_SUCCESS;
}

int main (void)
{
    double r = 0.0;
    double v = 0.0;
    double w = 2.4;
    double E1 = -2.3;
    double E2 = 0;
    double a = 0.7;
    gsl_odeiv2_system sys = {func, jac, 2, &r, &v, &w, &E1, &E2, &a};

    gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_x_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0);
    int i;
    double t = 0.0, t1 = 10000;
    double x[2] = {0.0, 0.0};

    for (i = 1 ; i<=10000; i++)
        {
            double ti = i*t1/10000;
            int status = gsl_odeiv2_driver_apply (d, &t, ti, x);

            if (status != GSL_SUCCESS)
                {
                    printf("error, return value%d\n", status);
                    break;
                }
            printf("%.5e %.5e %.5e\n", t, x[0], x[1]);
        }

    gsl_odeiv2_driver_free (d);
    return 0;
}

Answer 1

params 参数是指向某个任意数据结构的指针(地址/内存位置)。在 GSL 文档的示例中，他们的方程仅包含一个参数，这意味着可以只传递 double 字的地址。

但是，对于您的问题，您需要访问 6 个不同的参数。您无法使用相同的地址访问每个参数!

/* this doesn't work! */
double r  = *(double *)params;
double v  = *(double *)params;
double w  = *(double *)params;
double E1 = *(double *)params;
double E2 = *(double *)params;
double a  = *(double *)params;

由于所有地址都相同，因此您引用的是相同的号码。为了解决这个问题，您可以:将所有参数存储在长度为 6 的数组中，或者将它们存储在预定义的数据结构中。后一种方法更具可读性，因此我将对此进行演示。

首先定义一个数据类型来指定要存储的参数:

struct param_type {
    double r;
    double v;
    double w;
    double E1;
    double E2;
    double a;
};

现在，在 main 函数中创建这种类型的结构并存储参数的实际值:

struct param_type my_params = {r, v, w, E1, E2, a};

定义系统时，您存储指向该struct param_type的指针:

gsl_odeiv2_system sys = {func, jac, 2, &my_params};

要在 func 和 jac 中使用参数，只需从通用指针 (void *) 指向特定数据类型的指针 (struct param_type *):

struct param_type *my_params_pointer = params;

(请注意，在 C++ 中，必须使用显式强制转换来编写。)最后，您可以通过以下方式访问参数:

double r  = my_params_pointer->r;
double v  = my_params_pointer->v;
double w  = my_params_pointer->w;
double E1 = my_params_pointer->E1;
double E2 = my_params_pointer->E2;
double a  = my_params_pointer->a;

此处使用箭头 -> 而不是点 .，因为 my_params_pointer 是一个指针并且需要成为dereferenced使用前。