Search code examples
julianumerical-integration

Issues numerically integrating a multivariable function WRT a single variable in Julia (using hcubature)

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 enter image description here

Solution

  • 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.