Bivariate Gaussian Quadrature

A function implementing first, second and third order Gaussian qudrature in 2d cases with a test and illustration script attached.

For quadrature order over 3, symbolic toolbox is required to calculate the Gauss points and weights using Legendre Polynomials.

function quadrature = gaussianQuadrature2D(func, xmin, xmax, ymin, ymax, order)
%GAUSSIANQUADRATURE2D Perform a 2D Gaussian quadrature  over a specified range.
%   QUADRATURE = GAUSSIANQUADRATURE2D(FUNC, XMIN, XMAX, YMIN, YMAX, ORDER) 
% performs a 2D Gaussian quadrature numerical integration on a bivariate 
% function FUNC of order ORDER, over the range [XMIN, XMAX] and [YMIN, YMAX].
%   Algebra accuracy of ORDER-k Gaussian quadrature is (2*k - 1).
%
%   Inputs:
%   FUNC - a function handle to the integrand, supposed to be bivariate
%   XMIN, XMAX - the lower and upper bounds of the x range
%   YMIN, YMAX - the lower and upper bounds of the y range
%   ORDER - the order of Gaussian quadrature, determining the number of 
%   sample points and weights
%
%   Outputs:
%   QUADRATURE - the result of the numerical integration

A simple demo and test script:

% test_gauss_quad.m
clear
clc

% example 1
n = 10;
[xmin, xmax, ymin, ymax] = deal(-1.5, 2, -1, 0.5);
result = zeros(n, 5);
for k = 1:n
    func = @(x, y) (x + 0.5).^k .* (y - 0.3).^k + 1;
    
    % Gaussian quadrature with 1st/2nd/5th order
    gaussQuad1 = gaussianQuadrature2D(func, xmin, xmax, ymin, ymax, 1);
    gaussQuad2 = gaussianQuadrature2D(func, xmin, xmax, ymin, ymax, 2);
    gaussQuad3 = gaussianQuadrature2D(func, xmin, xmax, ymin, ymax, 5);
    
    % reference solution using built-in `integral2` function
    normalQuad = integral2(func, xmin, xmax, ymin, ymax);
    
    % logging result
    result(k, :) = [k, gaussQuad1, gaussQuad2, gaussQuad3, normalQuad];
end

% show result:
%  k-th order gaussian handles (2k-1) or lower polynomials accurately.
result = array2table(result); ...
result.Properties.VariableNames = ...
    {'Poly.Order', 'Gauss-1', 'Gauss-2', 'Gauss-5', 'Ref.Sol.'};
disp(result);

result:

    Poly.Order    Gauss-1    Gauss-2    Gauss-5    Ref.Sol.
    __________    _______    _______    _______    ________

         1        3.0844      3.0844     3.0844     3.0844 
         2        6.1433      9.3231     9.3231     9.3231 
         3        4.8815     -1.5406    -1.5406    -1.5406 
         4         5.402      17.027     19.903     19.903 
         5        5.1873     -15.115    -27.349    -27.349 
         6        5.2759      40.493     83.539     83.539 
         7        5.2393     -55.736    -189.11    -189.11 
         8        5.2544      110.78      504.8      504.8 
         9        5.2482     -177.37    -1309.3    -1309.3 
        10        5.2507      321.27     3534.2     3536.7