Gauss quadrature in 3D with change of variables

Specific Domains Numerics
1

Dear all,

I want to solve an intgeral of the following form
\int_0^1\int_x^1(x+y) dydx
using a self implemented Gauss quadrature.

The integral boundaries must be changed to [-1,1], which is done in the following way:
image

This code does not give the correct answer (the code is verified using the package QuadGK):

function f_example(x,y)
    return x + y 
end

function custom_gaussquad_change_of_variables(f, n)
    nodes, weights = gauss_quadrature_nodes_2d(n)
    integral = 0
   
    for i in 1:size(weights)[1]
        integral += weights[i] * f_example(1/2*nodes[i, 1]+1/2, (1-nodes[i,1])/2*nodes[i,2] + (1+nodes[i,1])/2) * 1/2 * (1-nodes[i,1])/2
    end 
     
    return integral
end

When I only need to change the boundary of one integral variable I do not have any problems. Does anyone understand what I do wrong and how I can obtain the correct answer?
Thanks!

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