MATLAB: Gaussian Quadrature Error Help

gaussian quadratureintegrationMATLABnumerical integration

Hi there, I have an assignment which asks to do Gaussian Quadrature for a 4th degree polynomial and discuss the error from the actual solution.
I'm a little confused however, as using 4 gauss points should give an exact solution for polynomials less than or equal to the degree 2n-1 = 7. Since I have a 4th degree polynomial here, why am I getting an error of around 3%? If I increase the number of gauss points to 9, I get an error in the order of 10^-8. Am I missing something major about Gaussian Quadrature? Thanks.
clear all
%% Exact solution
% Define the function as an inline function in x and y 
% Create an inline function
ff = @(x,y) 2.*x.^4 - x.^3 + 3.*y.^3 +y.^2 - 2.*x.*y + 5;
% Calculate the exact solution
ExactSol = integral2(ff,-1,1,-1,1);
% Plot the function
ezsurf(ff,[-1 1],[-1 1])
%% Gaussian Quadrature
n = 4;
xi = zeros(n,1);
eta = zeros(n,1);
weights = zeros(n,1);
evaluated = zeros(n,1);
multiplied = zeros(n,1);
% Location of the Gauss points
xi(1) = -1/(sqrt(3));
xi(2) = 1/(sqrt(3));
xi(3) = -1/(sqrt(3));
xi(4) = 1/(sqrt(3));
eta(1) = -1/(sqrt(3));
eta(2) = -1/(sqrt(3));
eta(3) = 1/(sqrt(3));
eta(4) = 1/(sqrt(3));
% Gauss weights
for i = 1:n
    weights(i) = 1;
end
% Evaulate the function at Gauss points, multiply by weights,
% then sum.
for j = 1:n
    evaluated(j) = 2*xi(j)^4 - xi(j)^3 + 3*eta(j)^3 +eta(j)^2 - 2*xi(j)*eta(j) + 5;
    multiplied(j) = evaluated(j)*weights(j);
end
GaussInt = sum(multiplied);
% Calculate the percentage error between the Gauss solution
% and the exact solution
Error = (GaussInt - ExactSol)/ExactSol*100;

Best Answer

I think you misiunderstand Gaussian quadrature.
You claim to have 4 points. WRONG. You have TWO points in x, and TWO in y, the result being a tensor product rule, formed using Gauss-Legendre nodes.
https://en.wikipedia.org/wiki/Gauss–Legendre_quadrature
An N point gaussian rule will be exact for polynomials of dergree 2*N-1. But your polynomial has degree 4. So why would you expect it to be exact? 2 points in each dimension means it will be exact only up to cubic functions in each variable.
If you increase the number of points to 9, thus a 3x3 tensor product rule, it should be exact (to within floating point trash) which just means you probably did something wrong there.
ff = @(x,y) 2.*x.^4 - x.^3 + 3.*y.^3 +y.^2 - 2.*x.*y + 5;
format long g
ffint = integral2(ff,-1,1,-1,1)
ffint = 
          22.9333333335287
That should be correct, to within the default tolerances integral2 returns. Just for kicks...
syms X Y
vpa(int(int(ff(X,Y),X,[-1,1]),Y,[-1,1]))
ans = 
22.933333333333333333333333333333
So I could probably have cranked down on integral2 just a bit more.
ffint = integral2(ff,-1,1,-1,1,'reltol',1e-15)
ffint = 
          22.9333333333334
I'm happpy with that, good enough, at least for government work. Let me now build 2 and 3 point tensor product rules.
leg2 = [-1 1]/sqrt(3);
w2 = [1 1];
leg3 = [-1 0 1]*sqrt(3/5);
w3 = [5 8 5]/9;
[xleg2,yleg2] = meshgrid(leg2);
[xleg3,yleg3] = meshgrid(leg3);
% perform the integration
ffgauss2 = w2*ff(xleg2,yleg2)*w2.'
ffgauss2 = 
          22.2222222222222

ffgauss3 = w3*ff(xleg3,yleg3)*w3.'
ffgauss3 = 
          22.9333333333333
So, as expected, the 3x3 node tensor product rule is exact. The 2x2 node rule is not. What a surprise. You should think about what I mean by a 3x3 rule as only 3 points in EACH dimension, NOT a 9 point rule. There is a HUGE difference.
I would also strongly suggest you learn to use MATLAB with vectors and arrays, perhaps as I did. Your code will improve by leaps and bounds when you start to do that.
As to what you did wrong? That part is up to you to diagnose, especially since I don't see where you wrote the code for the 3x3 rule.
Related Question