Question about 2D integral with a complex function f

[nhavt]

Dear @stevengj,

I would like to ask you a question regarding a 2D integral. Suppose that I have a complex function f as follows:

note:
x and y are scalars
A1, A2, A3, A4, and A5 are symmetric matrices
a and b are scalars
g1, g2, g3, and g4  are functions with parameters x, or y, or x and y

function f(x,y,A1,A2,A3,A4,A5,a,b) 
    C = g1(x)*A1 + g2(x)*A2+ g3(x,y)*A3 + A4 + g4(x,y)*A5
    Ω = pinv(C)
    return a*tr(Ω) + b*sum(Ω)
end

Suppose that parameters A1-A5, a, b g1-g4 are known, I would like to compute the following intergal:

I(f) = \int_{a_1}^{b_1}\int_{a_2}^{b_2} f(x,y,A1,A2,A3,A4,A5,a,b)dxdy

So, my question here is that would it be possible to compute I(f) using your package HCubarture with the complex function f? If possible, could you please tell me how to specify the function f in your package? Thank you very much in advance!

[mzaffalon]

From the README:

"Because hcubature is written purely in Julia, the integrand f(x) can return any vector-like object (technically, any type supporting +, -, * real, and norm: a Banach space). You can integrate real, complex, and matrix-valued integrands, for example."

[nhavt]

Hi @mzaffalon, thanks for your comment. However, from the readme, it mentions: "f should be a function f(x) that takes an n-dimensional vector x and returns the integrand at x". In my case, the function f(x, y, A1, A2, A3, A4, A5, a, b) has multiple parameters. The first two parameters (i.e. x and y) are related to the dimension of integral (2D), while the rest parameters are needed to compute the function at (x,y). Therefore, I don't know how to specify my function in hcubature. Do you have any suggestions for me?

[mzaffalon]

I didn't read your question properly, sorry. Yes, you can integrate it: just define your integrand as a function of x and y only

function f(x,y) 
    C = g1(x)*A1 + g2(x)*A2+ g3(x,y)*A3 + A4 + g4(x,y)*A5
    Ω = pinv(C)
    return a*tr(Ω) + b*sum(Ω)
end

The other constant must be assigned before you call hdubature(f, v1, v2): don't do the assignment in the global scope (i.e. outside a function) because the performance will be terrible. So you could do the following:

function compute_integral()
  A1 = ...
  A2 = ...
  function f(x,y) 
      C = g1(x)*A1 + g2(x)*A2+ g3(x,y)*A3 + A4 + g4(x,y)*A5
      Ω = pinv(C)
      return a*tr(Ω) + b*sum(Ω)
  end
  hcubature(f, v1, v2)
end

[nhavt]

Hi @mzaffalon, thanks very much for your help! Indeed, your second proposal is exactly what I am looking for! So, I will close the post!

[stevengj]

Alternatively, you could just do:

hcubature((x,y) -> f(x, y, A1, A2, A3, A4, A5, a, b), v1, v2)

[nhavt]

Dear @stevengj,

yes, that sounds great too! Thank you very much!