# Gauss quadrature in 3D with change of variables

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:

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

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!

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
ξ, 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


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

0.5000000000000002


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.