Gaussian Quadrature Error Help

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;  

Answer 1

[John D'Errico]

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.

[Cathal White]

Thank you for this answer, it clears up a lot for me. I wasn't trying to claim I was correct with my theory, so apologies if this was offensive. I'm an engineering student who needs this for doing some FE analysis. The lecturer briefly skimmed over the theory. Thanks again.

[John D'Errico]

Oh, god, you were not offensive. (I can be gruff at times.)

The point is, you can think of a 2x2 tensor product rule as I wrote it as integrating in x using 2 nodes, then integrating the result, which is now a function only of y also using a 2 node rule. The result is still only exact for 2*2-1=3, so cubic polynomials in each dimension.

Your problem is now of course to explain what I said, in your own words to discuss the error. But I think you have a better handle on it now.