Числено интегриране на сложна функция с Cubature C пакет

Искам да използвам пакета Cubature C, за да изпълня многомерен интеграл на сложна функция. Опитах се да го направя по следния начин за много простата функция 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