将数组传递到函数指针（在c）中？

Question

我正在尝试使用 Runge Kutta 解决二阶 ODE。我有 RK 方法的大纲，并且 ODE 本身已建立，但我无法继续。我在一个数组中有 ODE-s，我尝试将这个数组传递给 RK，但是在 f(dy[i]) 部分它给出了“预期的 double* 但键入 double”错误。当我尝试将两个 ODE-s 分成两个函数并将它们传递给 RK 时，它只适用于第一个函数。传递 dy 的元素而不是函数指针不是一种选择，因为 RK 本身应该能够求解包含任意数量变量的 ODE。

如何正确添加数组并使用它？ 这是我对应的一段代码：

double* harmOsc(double* p, double t, double* y, double* dy, int n)
{
    double k = p[0];
    double m = p[1];
    double x = y[0];
    double v = y[1];
    dy[0] = v;
    dy[1] = -k/m*x;
    return(dy);
}


double* HarmOsc1(double* dy, double t)
{
    return(dy);
}


void RK4(
    double t,       //independent variable
    double dt,      //stepsize

    double *y,      //variables
    double *dy,     //derivatives
    int n,          //number of equations
    double* (*f)(double*, double))
{

    int j;
    double* l = (double*)calloc(4*n,sizeof(double));
    for(j=0;j<n;j++)
    {
        l[0*n+j] = dt*f(dy[j],t);
        l[1*n+j] = dt*f(dy[j]+0.5*l[0*n+j],t+0.5*dt);
        l[2*n+j] = dt*f(dy[j]+0.5*l[1*n+j],t+0.5*dt);
        l[3*n+j] = dt*f(dy[j]+0.5*l[2*n+j],t+0.5*dt);
        y[j] = y[j] + (l[0*n+j] + 2*l[1*n+j] + 2*l[2*n+j] + l[3*n+j])/6;
    }
}

int main(int argc, char *argv[])
{
    double* p = (double*)calloc(2,sizeof(double));
    p[0] = 15; p[1] = 140;
    double* y = (double*)calloc(2,sizeof(double));
    y[0] = 12.4; y[1] = 1.1;
    double t=0;

    double* dy = (double*)calloc(2,sizeof(double));
    dy = harmOsc(p,t,y,dy,2);
    dy = HarmOsc1(dy,t);

    RK4(t,0.05,p,y,dy,2,&HarmOsc1);
    printf("%f, %f\n",dy[1],dy[0]);
}

HarmOsc1 是我调用的函数，因此它具有所需数量的参数。

当然还有我收到的警告：

RK3.c:55:13: 错误：‘f’的参数 1 的类型不兼容

RK3.c:55:13：注意：预期为“double *”，但参数类型为“double”

RK3.c:56:6: 错误：‘f’的参数 1 的类型不兼容

RK3.c:56:6：注意：预期为“双 *”，但参数的类型为“双”

RK3.c:57:6: 错误：‘f’的参数 1 的类型不兼容

RK3.c:57:6：注意：预期为“双 *”，但参数的类型为“双”

RK3.c:58:6：错误：‘f’的参数 1 的类型不兼容

RK3.c:58:6：注意：预期为“double *”，但参数类型为“double”

Answer 1

根据 cmets 中提供的说明，我们已经确定您的函数 RK4() 应该使用 Runge-Kutta 方法来数值求解以下形式的多个 ODE 中的每一个的一个值

y'(t) = f(y(t), t)

带初值条件的形式

y(t0) = y0

其中t 是一个标量，每个y 都是标量值，f() 对于所有 ODE 都是相同的（或至少由相同的 C 函数表示）。

我们进一步确定函数参数具有以下含义：

t    equivalent to t0 above
dt   the distance from t to the point at which the solutions are to be computed
y    points to an array in which to return the solutions
dy   points to the initial function values (y0 above) for all the ODEs
n    the number of ODEs to solve
f    the function `f()` above

你有各种各样的问题，但我认为大多数问题最终都是从你声明RK4()的参数f开始的。特别注意，ODE 的上述形式要求f() 接受一个与每个y 的值相同维度的值和另一个与每个y 的参数具有相同维度/数量的值作为参数，并且返回与每个y 的值相同维度的值。但是我们已经确定 ys 是 scalar 函数，每个函数都有 一个标量 参数，因此 f 应该接受两个 doubles 并返回一个 double .这将使我们得到这个声明：

void RK4(
    double t,       //initial point
    double dt,      //delta
    double *y,      //result values
    double *dy,     //initial values
    int n,          //number of equations
    double (*f)(double, double)) {

您不需要f 在一次运行中计算多个 ODE 的结果，因为您的 RK4() 迭代 ODE 并为每个 ODE 分别调用 f()。

接下来，让我们看看你的变量l。您正在动态分配足够的空间来存储每个输入 ODE 的所有四个 RK 常量（而不是释放它），但这没有用。每组 RK 常数只使用一次，因此在继续下一个方程时不需要记住以前的常数。因此，您只需要四个常量的空间，并且由于这是一个固定数字，您不需要动态分配。您甚至没有真正从使用数组中获得任何好处。我只会使用在循环体中声明的四个标量变量。

    int j;
    // double* l = (double*)calloc(4*n,sizeof(double));
    for (j = 0; j < n; j++) {
        double k1 = dt * f(dy[j], t);
        double k2 = dt * f(dy[j] + 0.5 * k1, t + 0.5 * dt);
        double k3 = dt * f(dy[j] + 0.5 * k2, t + 0.5 * dt);
        double k4 = dt * f(dy[j] + k3, t + dt);  // <-- note corrections

请注意，现在声明 f 的参数和返回类型符合 RK 方程对输入 ODE 形式的要求。

最后，计算结果：

        y[j] = dy[j] + (k1 + 2 * k2 + 2 * k3 + k4) / 6;

请注意最后一行和您的版本之间的本质区别——RK 计算增量 y，但您的初始 y 值在 dy 中，而不是在 y 中。

就是这样：

    }
}

此时我观察到，我在解决您的问题时遇到的最大问题来自于不了解您尝试做的事情的细节。这源于几件事，其中很大

• 您的变量和参数名称很短且无法表达，更糟糕的是，就它们似乎的含义而言，它们中的一些实际上代表的东西与其名称所暗示的不同

• 代码文档很少，而且显然不正确。这似乎与前一点有点相符，因为鉴于有限（但不完全正确）的参数文档，函数参数命名更有意义。

• 代码本身的帮助有限，因为引发问题的正是缺陷。

带回家：编写好的代码文档。使用完整的句子。描述每个函数对每个参数的期望，以及它承诺用它们做什么。描述返回值。记录函数执行的任何错误处理。当你这样做的时候，试着像一个想要使用你的函数的人一样思考，但不能咨询它的实现来看看它做了什么——那个人需要知道什么？事实上，我建议在编写每个函数的实现之前为每个函数编写文档。请根据您的需要更新文档，但文档可以帮助您保持进度，并且首先编写它们可以帮助您编写，就好像您看不到实现一样。

Answer 2

函数f 期待double *，即指向double 的第一个参数的指针。但是，在您调用f 的每个地方，您传递的是double 值，而不是指向double 的指针，也不是数组（衰减为指向第一个值的指针）：

    // here ------------v
    l[0*n+j] = dt*f(  dy[j],                 t);
    l[1*n+j] = dt*f(  dy[j]+0.5*l[0*n+j],    t+0.5*dt);
    l[2*n+j] = dt*f(  dy[j]+0.5*l[1*n+j],    t+0.5*dt);

如果这个函数的第一个参数需要是一个数组，那么传递一个数组，而不是单个值。

Answer 3

f 期望第一个参数是 double*，而您使用 dy[j] 调用 f，它是一个 double。 只需将 dy[j] 更改为 dy + j 即可。

或

你可能需要一个临时变量

double tempv;

然后在循环内，

tempv = dy[j];
...
tempv = dy[j] + ...;

用tempv的地址调用f，

f(&tempv,...)

Answer 4

在 C 中，不建议返回数组。这是可能的，但管理负责释放此内存的代码很麻烦。此外，它违背了一些 C 哲学，即生成高性能代码不必要地重复生成本质上相同的数组，只是在此后不久将其销毁。在脚本语言中对算法进行原型设计时，您可以这样做。

因此最好让导函数有签名

void f( double * dy, double * y, double t);

其中维度和常量是外部常量（或将参数数组添加到参数列表中）

然后你在 rk4 循环中使用它

f(k1,  y, t      );
for(i=0; i<n; i++) yt[i] = y[i] + 0.5*h*k1[i];
f(k2, yt, t+0.5*h);
for(i=0; i<n; i++) yt[i] = y[i] + 0.5*h*k2[i];
f(k3, yt, t+0.5*h);
for(i=0; i<n; i++) yt[i] = y[i] +     h*k3[i];
f(k4, yt, t+h);
for(i=0; i<n; i++) y[i] = y[i] + h/6*(k1[i]+2*(k2[i]+k3[i])+k4[i]);

当然，您需要确保每次集成也只分配一次k1,k2,k3,k4,yt。这可以通过 3 种方式发生：

• 积分过程包含完整的循环，必要时返回一个点列表，
• 有些结构包含在第一个集成步骤之前构建的工作数组，
• 或者您使用更新的动态堆栈分配，它也不应该有一个禁止的分配惩罚。