2D Integrals : Gaussian Quadrilateral Example

Table of Contents

A first example on the reference quadrilateral

Let's consider the function defined on the reference quadrilateral with the vertices
We want to compute numerically the double integral
by using the 2D Gauss formula which is similar to the 1D case but it involves the 2D Gauss points with weights
. where N is the total number of Gauss points.

Compute the 2D Gauss points on the reference element

First we compute the appropriate Gauss points in the reference quadrilateral. They are obtained from the 1D Gauss points in [-1,1] in both axes x and y and are generated by the Matlab function gaussValues2DQuad which can be found at the end of this script.
Let Nax denote the number of 1D Gauss points in each of the axes x,y. Then the number of N of 2D Gauss points in is N=Nax^2. Notice that if is a polynomial function of degree n < 2·Nax+1 in each of the variables, then the 2D Gauss integration formula returns an exact value of the integral. Since in our case , to get an exact result it is sufficient to consider the case
Nax=2; %order of the Gaussian quadrature
fprintf('N = %d, w= %f %f %f %f\n', 2*Nax,w);
N = 4, w= 1.000000 1.000000 1.000000 1.000000
ptGaussRef % +/- 1/sqrt(3)
ptGaussRef = 4×2
-0.5774 -0.5774 -0.5774 0.5774 0.5774 -0.5774 0.5774 0.5774
% Display Gauss points in the reference quadrilateral [-1,1]x[-1,1]
vertices=[-1 -1; 1,-1; 1,1; -1,1; -1 -1];
plot(vertices(:,1), vertices(:,2));
hold on;
hold off;
F=@(x,y) (x.^2).*y.^2;
Integ = 0.4444

General case: A quadrilateral domain

Let's consider the function
and the quadrilateral with the vertices
We want to compute the following integral by using the formula of change of variables
where is the absolute value of the determinat of the Jacobian matrix of the change of coordinates.
Below, you have an slide showing the change of variables needed to relate the reference quadrilateral [-1,1]x[-1,1] with a general one. The are known as shape functions associated to each vertice and are explained below:
are functions associated to the respective vertex of the refrence quadrilateral and satisfy
if and


Now consider our function and we code in Matlab the above change of variable formulas.
% The present function
F =@(x,y) (x.^2)- y.^2;
% The domain vertices
% We use Matlab *implicit function* definition:
% Shape functions
% Shape function derivatives
dPsi11=@(xi,eta) -(1-eta)/4;
dPsi21=@(xi,eta) (1-eta)/4;
dPsi31=@(xi,eta) (1+eta)/4;
dPsi41=@(xi,eta) -(1+eta)/4;
dPsi12=@(xi,eta) -(1-xi)/4;
dPsi22=@(xi,eta) -(1+xi)/4;
dPsi32=@(xi,eta) (1+xi)/4;
dPsi42=@(xi,eta) (1-xi)/4;
% Gradient matrix
Jacb =@(xi,eta) [dPsi11(xi,eta), dPsi21(xi,eta),dPsi31(xi,eta),dPsi41(xi,eta);...
dPsi12(xi,eta), dPsi22(xi,eta),dPsi32(xi,eta),dPsi42(xi,eta)];

Compute the corresponding Gaussian points on the domain

As we know from above, the change of variables from the reference quadrilateral to a general one is:
Screenshot 2023-11-23 at 11.24.22.png
evaluate Shape functions on Gaussian reference points
xi = ptGaussRef(:,1);
eta = ptGaussRef(:,2);
evalPsi1 = Psi1(xi,eta);
evalPsi2 = Psi2(xi,eta);
evalPsi3 = Psi3(xi,eta);
evalPsi4 = Psi4(xi,eta);
% change these points to the new quadrilateral(from the change of variables function)
ptGaussDomain = evalPsi1*v1+evalPsi2*v2+evalPsi3*v3+evalPsi4*v4;
% Draw Gauss points in the present domain
plot(vertices(:,1), vertices(:,2));
hold on;
hold off;

Compute the Jacobian terms

vertex matrix
v = [v1;v2;v3;v4];
% evaluate Jacobian contribution for each Gauss point
for i=1:size(xi,1)
evalDetJacb(i) = abs(det(Jacb(xi(i),eta(i))*v));

Compute the integral value according Gauss formula

%evaluate the function on the domain points
% Finally, apply Gauss integration formula
for i=1:size(ptGaussDomain,1)
integ = suma
integ = 60.0000
Warning 1. As we know from the theory, the area of any elementary domain in the plane (x,y) and, in particular, of a quadrilateral , can be calculated via the double integral
Note that if we apply the above MATLAB code to this case (or the case of any constant function) by writing the subintegral function as
F=@(x,y) 1
F=@(x,y) 1.
the program will return an error message: "Index exceeds the number of array elements".
One of the simplest solutions would be to write instead, for example
F=@(x,y) ones(size(x));
F=@(x,y) x./x
but then we must ensure that the function F will not be evaluated at points with zero x-coordinate.
The definitive solution is
F=@(x,y) (x.^2+y.^2+1)/(x.^2+y.^2+1)
Warning 2. In the above list v=[v1; v2; v3; v4] the vertices v1,...,v4 of the quadrilateral must be arranged in the counterclockwise or clockwise order around . Namely, for our choice
v1=[0,0]; v2=[5,-1]; v3=[4,5]; v4=[1,4];
we can write
v=[v1; v2; v3; v4] or v=[v4; v3; v2; v1] or v=[v2; v3; v4, v1]
(this will give the same value of the integral),
but we cannot write
v=[v1; v3; v2; v4] or v=[v4; v2; v1; v3],
as it will lead to a wrong value of the determinant evalDetJacb(i) of the Jacobian matrix.

Exercise 1: Build the integQuad function

From the previous example, one can build a function that returns the value of the integral of an inline function defined on a quadrilateral domain (see at the end of this document).
function integralValue=integQuad(F,vertices,N)
(Hint: find the answer at the end of this document)

Application: Integration over a mesh

Let's consider the integral of
in a general rectangular domain , in particular we know that
Instead of integrating only in one quadrilateral, we want to show how to approximate the integral value using a mesh subdivision of the domain
The final integral will be the sum of the integrals in each small quadrilateral. Therefore, we increase the precission on the final value.
F=@(x,y) 2*exp(x+y.^2).*y; %present function
xmin=0; % rectangle dimensions
% -----------------------------------------------
% In order to create a rectangular mesh in Matlab
% we can make the following procedure
% -----------------------------------------------
numDiv=4; %number of subdivisions
xi=xmin:hx:xmax; %all points in the x-axis
eta=ymin:hy:ymax; %all points in the y-axis
[X,Y] = meshgrid(xi,eta); % create a grid of points
Z = F(X,Y); % function values
surf(X,Y,Z); % optional: is just to visualize the result
title('Plot of F(x,y) in the mesh domain')
Hint. Here we plot the surface as the approximation of the continous function using the mesh
[elem,vertex] = surf2patch(X,Y,Z); % the main variables
elem = 16×4
1 6 7 2 2 7 8 3 3 8 9 4 4 9 10 5 6 11 12 7 7 12 13 8 8 13 14 9 9 14 15 10 11 16 17 12 12 17 18 13
vertex = 25×3
0 0 0 0 0.2500 0.5322 0 0.5000 1.2840 0 0.7500 2.6326 0 1.0000 5.4366 0.2500 0 0 0.2500 0.2500 0.6834 0.2500 0.5000 1.6487 0.2500 0.7500 3.3803 0.2500 1.0000 6.9807
numElem=size(elem,1); %total number of elements
numVert= size(vertex,1); % total number of vertices
fprintf('num Elements =%d, num vertices = %d \n',numElem,numVert)
num Elements =16, num vertices = 25


Fig. Visualization of the numbering for vertex and elements for a 4x4 mesh.

Structure of the Mesh

In general a mesh is defined by two tables: vertex and elem.
vertex: is a matrix with 3 columns and numVert rows, where numVert is the number of grid points.
For each row:
elem: is a numElem×4 matrix, where numElem is the number of elements For each row:
N=2; %number of Gauss points = NxN
for i=1:numElem %compute the integral on each element
n1=elem(i,1); %first vertex of the i-th element
n2=elem(i,2); %second vertex of the i-th element
n3=elem(i,3); %third vertex of the i-th element
n4=elem(i,4); %fourth vertex of the i-th element
actualIntegVal= (exp(1)-1)^2 %exact known value
actualIntegVal = 2.9525
errorInt=abs(actualIntegVal-integralTot) %absolute error
errorInt = 2.9420e-04

Additional Functions

Gauss points in 2D

They are the same as in the 1D version, but now using the points in both x and y axis.
function [w2D,pt2D]=gaussValues2DQuad(n)
switch (n)
case 1
w=2; pt=0;
case 2
w=[1,1]; pt=[-1/sqrt(3), 1/sqrt(3)];
case 3
w=[5/9, 8/9, 5/9]; pt=[-sqrt(3/5), 0, sqrt(3/5)];
case 4
w=[(18+sqrt(30))/36, (18+sqrt(30))/36, (18-sqrt(30))/36, (18-sqrt(30))/36];
pt=[-sqrt(3/7-2/7*(sqrt(6/5))), sqrt(3/7-2/7*(sqrt(6/5))),...
-sqrt(3/7+2/7*(sqrt(6/5))), sqrt(3/7+2/7*(sqrt(6/5)))];
case 5
w=[(322+13*sqrt(70))/900, (322+13*sqrt(70))/900, 128/225,...
pt=[-1/3*sqrt(5-2*sqrt(10/7)), 1/3*sqrt(5-2*sqrt(10/7)),0,...
-1/3*sqrt(5+2*sqrt(10/7)), 1/3*sqrt(5+2*sqrt(10/7))];
error('No data are defined for this value');
for i=1:n
for j=1:n
pt2D=[pt2D; [pt(i),pt(j)]];

Quadrilateral integration function: integQuad

Copy this function in a new file named integQuad.m
function valInteg = integQuad(F,vertices,N)
% Shape functions
% Shape function derivatives
dPsi11=@(x,y) -(1-y)/4;
dPsi21=@(x,y) (1-y)/4;
dPsi31=@(x,y) (1+y)/4;
dPsi41=@(x,y) -(1+y)/4;
dPsi12=@(x,y) -(1-x)/4;
dPsi22=@(x,y) -(1+x)/4;
dPsi32=@(x,y) (1+x)/4;
dPsi42=@(x,y) (1-x)/4;
% Gradient matrix
Jacb =@(x,y) [dPsi11(x,y), dPsi21(x,y),dPsi31(x,y),dPsi41(x,y);...
dPsi12(x,y), dPsi22(x,y),dPsi32(x,y),dPsi42(x,y)];
% evaluate Shape functions on Gaussian reference points
xi = ptGaussRef(:,1);
eta = ptGaussRef(:,2);
evalPsi1 = Psi1(xi,eta);
evalPsi2 = Psi2(xi,eta);
evalPsi3 = Psi3(xi,eta);
evalPsi4 = Psi4(xi,eta);
% from the change of variables function
ptGaussDomain = [evalPsi1,evalPsi2,evalPsi3,evalPsi4]*vertices;
% evaluate Jacobian contribution for each point
for i=1:size(xi,1)
evalDetJacb(i) = abs(det(Jacb(xi(i),eta(i))*vertices));
%evaluate the function on the domain points
% Finally, apply Gauss formula
for i=1:size(ptGaussDomain,1)
valInteg = suma;

Exercise 2:

Increase the number of subdivisions and see how the error evolves.
(c) Numerical Factory 2022