# Gauss quadrature in 3D with change of variables
**URL:** https://discourse.julialang.org/t/gauss-quadrature-in-3d-with-change-of-variables/113486
**Category:** Numerics
**Created:** [April 25, 2024, 7:52am UTC](https://discourse.julialang.org/t/gauss-quadrature-in-3d-with-change-of-variables/113486 "2024-04-25T07:52:08Z")
**Posts on this page:** 2
**Page:** 1
### Author:  [@Lisette\_de\_Bruin](https://discourse.julialang.org/u/Lisette_de_Bruin)
#### Post date: [April 25, 2024, 7:52am UTC](https://discourse.julialang.org/t/gauss-quadrature-in-3d-with-change-of-variables/113486/1 "2024-04-25T07:52:08Z")
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):
```julia
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!
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [April 25, 2024, 12:45pm UTC](https://discourse.julialang.org/t/gauss-quadrature-in-3d-with-change-of-variables/113486/2 "2024-04-25T12:45:43Z")
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:
```julia
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
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](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-changing-the-type-of-a-variable). If you know you want a `Float64` result, you can use `integral = 0.0` to initialize it.