Gauss quadrature in 3D with change of variables

Specific Domains Numerics
2

Your basic idea is correct, but you have a bug/typo in your code somewhere, and you also didn’t supply gauss_quadrature_nodes_2d so it is not runnable. I found your code a bit hard to read so I re-implemented it, and mine seems to work fine:

using QuadGK
function custom_gaussquad(f, n)
    ξ, w = QuadGK.gauss(n, -1, 1) # 1d gauss rule of order n
    integral = 0.0 # not 0, for type stability
    for (ξ₁,w₁) in zip(ξ, w), (ξ₂,w₂) in zip(ξ, w) # tensor product of 1d rules
        x = (ξ₁+1)/2
        y = ((1-x)*ξ₂ + (1+x))/2
        integral += (w₁ * w₂) * f(x, y) * (1-x)/4
    end
    return integral
end

which gives the correct answer:

julia> custom_gaussquad((x,y) -> 1, 11) # area of triangle
0.49999999999999983

julia> custom_gaussquad((x,y) -> x+y, 11)
0.5000000000000002

julia> quadgk(x -> quadgk(y -> x+y, x, 1)[1], 0,1)[1]
0.49999999999999994

julia> custom_gaussquad((x,y) -> exp(x^2 + sin(x*y)), 11)
0.8215954352465259

julia> quadgk(x -> quadgk(y -> exp(x^2 + sin(x*y)), x, 1)[1], 0,1)[1]
0.8215954352465216

PS. Note that integral = 0 is type unstable. If you know you want a Float64 result, you can use integral = 0.0 to initialize it.

1 Like