Численное интегрирование сложной функции с помощью пакета Cubature C

Я хочу использовать пакет C Cubature для выполнения многомерного интеграла сложная функция. Я попытался сделать это следующим образом для очень простой функции f (x, y) = x + y * i над квадратом [0,1] x [0,1]. Точный результат - 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?


c math complex-numbers numerical-methods
person Šatov    schedule 28.05.2013    source источник
comment
Я ничего не знаю об этом пакете, но мне непонятно, как fval[0] = x[0]+x[1]*I; мог успешно упаковать комплексное число в один double   -  person AakashM    schedule 28.05.2013
comment
тег R на самом деле? В R (язык) есть кубатурный пакет.   -  person baptiste    schedule 28.05.2013
comment
Да, есть! Вы можете найти его здесь   -  person Šatov    schedule 28.05.2013
comment
@AakashM Вы абсолютно правы.   -  person Šatov    schedule 28.05.2013
comment
@RiccardoCavallari Я знаю, но я хотел сказать следующее: относится ли этот вопрос к пакету R или нет?   -  person baptiste    schedule 28.05.2013
comment
@baptiste Это не строго связано с R, поэтому я удалил тег, если это было вашей точкой зрения, но я подумал, что люди, использующие Cubature в R, могли столкнуться с той же проблемой.   -  person Šatov    schedule 28.05.2013
comment
Я не испытываю сильных чувств к вашему тегу R, но вы могли бы более четко объяснить, как он связан с вопросом.   -  person baptiste    schedule 28.05.2013
comment
Почему бы просто не объединить реальную и мнимую части по отдельности? У него не будет оптимальной производительности, потому что некоторые вычисления будут выполняться дважды, но это будет не медленнее, чем вдвое медленнее, чем вы могли бы сделать.   -  person sigfpe    schedule 29.05.2013

Ответы (2)


arrow_upward
3
arrow_downward

2 вопроса:

1) Вы объявляете double val, но хотите, чтобы это было double val[2]. Ваш исходный val, являющийся double в cimag (val), всегда будет возвращать 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" для тестирования прямо сейчас.

person chux - Reinstate Monica    schedule 05.06.2013

arrow_upward
0
arrow_downward

Когда вы определяете val, сделайте следующее:

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

где integrand_fdim равно unsigned fdim - первый аргумент вашего подынтегрального выражения f.

P.S .: ответ @chux производит Bus 10 или segfault ...

person MsTais    schedule 23.03.2018