2

我想使用Cubature C 包来执行复杂函数的多维积分。对于正方形 [0,1]x[0,1] 上的非常简单的函数 f(x,y) = x + y*i，我尝试按以下方式执行此操作。确切的结果是 0.5 + 0.5i。

#include <stdio.h>
#include <math.h>
#include <complex.h>
#include "../cubature.h"

int f(unsigned ndim, const double *x, void *fdata, unsigned fdim, double *fval);

int main(void)
{
    double xmin[2] = {0,0}, xmax[2] = {1,1}, val, err;
    hcubature(1, f, NULL, 2, xmin, xmax, 0, 0, 1e-3, ERROR_PAIRED, &val, &err);
    printf("Computed integral = %f+%fi +/- %f\n", creal(val),cimag(val), err);
}

int f(unsigned ndim, const double *x, void *fdata, unsigned fdim, double *fval) 
{
    fval[0] = x[0]+x[1]*I;
    return 0;
}

但我得到的是Computed integral = 0.500000+0.000000i +/- 0.000000，那是虚部没有出现。如果我放一个纯虚被积函数（例如 x*i），我总是得到 0 作为结果。

我究竟做错了什么？你知道在 C 中计算复杂函数积分的更好方法吗？

4

2 回答 2

3

2个问题：

1）您正在声明double val，但您希望它是double val[2]。您的原始 val 是doublecimag(val) 中的 a 将始终返回 0.0。改成

double xmin[2] = {0,0}, xmax[2] = {1,1}, val[2] /* not val */, err;
hcubature(2 /* not 1 */, f, NULL, 2, xmin, xmax, 0, 0, 1e-3, ERROR_PAIRED, /* & */val, &err);
printf("Computed integral = %f+%fi +/- %f\n", val[0], val[1], err);

2) 看来您在 中的计算不f(x,y) = x + y*i正确f()。 fval[0]是双重的，无法容纳您尝试分配的复杂答案。

int f(unsigned ndim, const double *x, void *fdata, unsigned fdim, double *fval) {
    fval[0] = x[0];
    fval[1] = x[1];
    return 0;
}

抱歉，我现在没有要测试的“../cubature.h”。

于 2013-06-05T05:08:53.583 回答
0

定义 val 时，请执行以下操作：

double *val=(double *) malloc(sizeof(double) * integrand_fdim);

其中integrand_fdim等于unsigned fdim- 你的被积函数的第一个参数f。

PS：@chux 的回答会产生 Bus 10 或 segfault ......

于 2018-03-23T19:27:35.313 回答