# MATLAB: 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 functionff = @(x,y) 2.*x.^4 - x.^3 + 3.*y.^3 +y.^2 - 2.*x.*y + 5;% Calculate the exact solutionExactSol = integral2(ff,-1,1,-1,1);% Plot the functionezsurf(ff,[-1 1],[-1 1])%% Gaussian Quadraturen = 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 pointsxi(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 weightsfor 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);endGaussInt = sum(multiplied);% Calculate the percentage error between the Gauss solution% and the exact solutionError = (GaussInt - ExactSol)/ExactSol*100;  ``

• 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.
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 gffint = 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 Yvpa(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 integrationffgauss2 = 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.