hcubature
keeps returning zeros for me.
I have a multivariable function that I would like to integrate with respect to one variable. As a test, I use quadgk
, which can only be used for 1 variable (hardcoding the other variables), and it works.
However, when using hcubature
my integration returns zeros.
here is my minimal non-working example
psi_1D_PIB(x, n, L) = sqrt(2/L) * sin(n * π * x / L)
psi_conj_PIB(x, n, L) = conj(psi_1D_PIB(x, n, L)) * psi_1D_PIB(x, n, L)
psi_conj_PIB(v) = psi_conj_PIB(v...) # as per hcubature documentation, send an array
println(hcubature(psi_conj_PIB, (1.,1, 5.), (5.,1,5.))) # integrate x from 1 to 5, n from 1 to 1, and L from 5 to 5.
# first tuple is LH bounds, second is RH integration bounds
output is 0.0, 0.0)
The missing left hand parenthesis does not appear to be a typo
Using guadgk
to calculate this essentially single variable integral,
psi_1D_PIB(x) = sqrt(2/5) * sin(1 * π * x / 5)
psi_conj_PIB(x) = conj(psi_1D_PIB(x)) * psi_1D_PIB(x)
println(quadgk(psi_conj_PIB, 1., 5.))
output is (0.9513653457281316, 2.5028679129235343e-10)
For additional information, this is the function being integrated
You must use:
julia> println(hcubature(x->psi_conj_PIB(first(x),1,5),[1.],[5.]))
(0.9513653457281317, 2.5028645822544604e-10)
You have a function f(x,y,z) and want to compute its integral over the interval x in [1,5] for given values y=1, z=5.
The right way to do that is to restrict your R^3 function f(x,y,z) to a R function: x->g(x) = f(x,1,5) and to compute the integral of g over the [1,5] interval. In Julia this is written
x->psi_conj_PIB(first(x),1,5)
Note: here you must use an extra first(x)
to get a scalar from the x vector of dimension one, because your R^n domain of integration is in fact a simple R domain defined by lower bound = [1], upper bound = [5].
Note 2: this is not equivalent to compute an integral over the region (x,y,z) in [1,5]x[1,1]x[5,5]. This R^3 region has a null "volume": (5-1)x(1-1)x(5-5) = 4x0x0 = 0, thus whatever the function you want to integrate its integral over this region is null. This is what you experimented in your example.