Gaussian Quadrature Error Help
20 views (last 30 days)
Show older comments
Cathal White
on 19 Nov 2020
Commented: John D'Errico
on 19 Nov 2020
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;
0 Comments
Accepted Answer
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)
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]))
So I could probably have cranked down on integral2 just a bit more.
ffint = integral2(ff,-1,1,-1,1,'reltol',1e-15)
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.'
ffgauss3 = w3*ff(xleg3,yleg3)*w3.'
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
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.
More Answers (0)
See Also
Categories
Find more on Polynomials in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!