4.4

4.4 | 10 ratings Rate this file 30 Downloads (last 30 days) File Size: 1.84 KB File ID: #9230

Gaussian Quadrature for Triangles

by

 

03 Dec 2005 (Updated )

Compute Gauss nodes and weights for a triangle

| Watch this File

File Information
Description

This script computes Gauss nodes and weights for numerically integrating a fuction over an arbitrary triangular domain. Note: This method uses the collapsed square rather than the general cubature case.

MATLAB release MATLAB 7 (R14)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (16)
14 Jun 2012 Cássio

Hey Greg, great code!
Could you tell me a book/paper to read more about?
Thanks

01 Apr 2012 sfilop  
28 Jan 2011 Grzegorz Knor

Ok, it was my mistake, should be:
f=@(x,y) x.^2.*exp(y);
instead of:
f=@(x,y) x^.2.*exp(y);

I wrote a function (slower) from which you can compare the results:
================================
function sol = int_tri(fun,p)
xa = min(p(:,1));
xb = median(p(:,1));
xc = max(p(:,1));
ya = p(p(:,1)==xa,2);
yb = p(p(:,1)==xb,2);
yc = p(p(:,1)==xc,2);
if length(unique(p(:,1)))==3 % two triangles
ys = (ya-yc)/(xa-xc)*xb + (xa*yc-xc*ya)/(xa-xc);
y = [([xa 1; xb 1]\[ya; min(yb,ys)])' % bottom line 1
([xa 1; xb 1]\[ya; max(yb,ys)])' % top line 1
([xb 1; xc 1]\[min(yb,ys); yc])' % bottom line 2
([xb 1; xc 1]\[max(yb,ys); yc])'];% top line 2
y_lo1 = @(x)y(1,1)*x+y(1,2);
y_lo2 = @(x)y(3,1)*x+y(3,2);
y_hi1 = @(x)y(2,1)*x+y(2,2);
y_hi2 = @(x)y(4,1)*x+y(4,2);
fh1 = @(x) quadl(@(y)fun(x,y),y_lo1(x), y_hi1(x));
fh2 = @(x) quadl(@(y)fun(x,y),y_lo2(x), y_hi2(x));
sol = quadl(@(X)arrayfun(fh1,X), xa,xb) + quadl(@(X)arrayfun(fh2,X), xb,xc);
else % one triangle
y = [([xa 1;xc 1]\[min(ya); min(yc)])'
([xa 1;xc 1]\[max(ya); max(yc)])'];
y_lo1 = @(x)y(1,1)*x+y(1,2);
y_hi1 = @(x)y(2,1)*x+y(2,2);
fh1 = @(x) quadl(@(y)fun(x,y),y_lo1(x), y_hi1(x));
sol = quadl(@(X)arrayfun(fh1,X), xa,xc);
end

25 Nov 2010 Grzegorz Knor

Why my results depend on N?
For example:
f=@(x,y) x^.2.*exp(y);

[X,Y,Wx,Wy]=triquad(1,[0 0; 0 2; 2 0]);
Q=Wx'*feval(f,X,Y)*Wy

[X,Y,Wx,Wy]=triquad(10,[0 0; 0 2; 2 0]);
Q=Wx'*feval(f,X,Y)*Wy

[X,Y,Wx,Wy]=triquad(20,[0 0; 0 2; 2 0]);
Q=Wx'*feval(f,X,Y)*Wy
?

10 Feb 2010 saurabh gupta  
18 Aug 2008 DAR .

Seems good. Too much simple thought.

25 May 2008 meysam mousavi  
04 Mar 2008 Im Serious

Does not integrate correctly

14 Aug 2007 YALA Hakim

I need this file for my code of FEM

28 Mar 2007 zhenzhen Jin  
08 Nov 2006 Greg von Winckel

First, the area would not be 1, but rather 1/2. Second, you need to define a function like f=@(x,y) 1+0*(x.*y) to make sure it returns something with the correct dimension. Alternately, you could just sum the weights to get the area.

30 Oct 2006 marki bugs

say vertices are (0,0), (1,0), (0,1)
I see this function & example works well, but besides f=(x,y) 1. I understand if I insert f=(x,y) 1, i should get the area of the triangle back. But I don't, am I not doing something wrong?

12 Aug 2006 John D'Errico

Garrett defined his function improperly. The lack of a dot can be important in matlab.

f=@(x,y) x.^2 + exp(y);

09 Aug 2006 Garrett Baird

I must be doing something wrong, because I can't get this to work for me. For example

[X,Y,Wx,Wy]=triquad(8,[0 0; 0 2; 2 0]);
f=@(x,y) x^2 + exp(y);
Q=Wx'*feval(f,X,Y)*Wy;

returns Q = 7.63, but the real answer is exp(3) - 5/3 which is around 5.7224. Increasing the number of quadrature points just makes Q higher.

06 Dec 2005 Greg von Winckel

Thanks for catching that John. The comments have been corrected.

05 Dec 2005 John D'Errico

A nice code, but with an error in the documentation -
the example usage fails. I'll assume that the author really
intended to write:

Q=Wx'*feval(f,X,Y)*Wy;

A quick test verifies that this did give the exact integral
on a simple linear function over a simple triangle.

I'll still rate it a 5, and hope the author chooses to fix the
documentation.

Updates
07 Dec 2005

Missing transpose symbol in comments.

Contact us