Gaussian Quadrature Error Help

20 views (last 30 days)
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;

Accepted Answer

John D'Errico
John D'Errico on 19 Nov 2020
Edited: John D'Errico on 19 Nov 2020
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 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.
  2 Comments
Cathal White
Cathal White on 19 Nov 2020
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
John D'Errico on 19 Nov 2020
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Polynomials in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!