A question regarding nodes and weights in hcubature

Dear @stevengj

I would like to ask you a question when computing ∫∫g(x)*g(y)f(x,y)dxdy. I hope it is fine to post my question here. Below is my problem.

using QuadGK, HCubature
# get a full grid
function fullGrid(x1,x2,w1,w2)
    pW = vec(collect.(Iterators.product(w1,w2)))
    pX = vec(collect.(Iterators.product(x1,x2)))
    lp = length(pW)
    w = zeros(lp)
    x = zeros(lp,2)
    for i = 1:lp
        w[i] = pW[i][1]*pW[i][2]  
        x[i,1] = pX[i][1]
        x[i,2] = pX[i][2]       
    end
    return x,w
end
g(x) = 1/(1+x^2) # this is my weight function
h(x) = 1 # I simplify the function h(x). In my case, h(x) is quite complex and is very costly to compute.

I would like to compute ∫∫g(x)*g(y)f(x,y)dxdy in a hypercube [a, b]^2 where a,b≥0. Here I use two approaches based on your two packages (QuadGK and HCubature).

function test(n,lb,ub)
    # Approach 1: 
    # get quadrature nodes for each dimension with a weight function g(x)
    x1, w1 = QuadGK.gauss(x ->g(x), n, lb[1], ub[1], rtol=1e-10)
    x2, w2 = QuadGK.gauss(x ->g(x), n, lb[2], ub[2], rtol=1e-10)
    # obtain a full grid
    x, w = fullGrid(x1,x2,w1,w2)
    s1 = sum(h(x) .* w) # I made a mistake here. f ---> h
    f(x) = g(x[1])*g(x[2])*h(x)
    # Approach 2: using HCubature as a benchmark
    s2 = hcubature(x -> f(x), lb, ub)
return s1, s2[1]
end

Below are the results with the two approaches:

n = 20 # the number of sample points
lb1 = [0,0]
ub1 = [1,1]
# case 1
res1 = test(n,lb1,ub1)
# res1 = (0.6166469119120699, 0.616850275222724)

In this case, it looks fine with Approach 1.

# case 2:
lb2 = [6,6]
ub2 = [10,10]
res2 = test(n,lb2,ub2)
# res2 =(0.004287633663977574, 0.004287633663975462)

Approach 1 still works in this case when changing the intervals even though the weight function g(x) is a rapid decrease when x is large. Since I would like to get nodes and weights to compute∫∫g(x)*g(y)f(x,y)dxdy, could I extract nodes and weights with hcubature given the weight functions g(x) and g(y), and the function f(x,y) is very costly to compute?

Thank you very much for your help. I am sorry if my question is very basic since I am a beginner at numerical integration.

JuliaMath / HCubature.jl

A question regarding nodes and weights in hcubature #38