Code covered by the BSD License

### Highlights from Gaussian Quadrature for Triangles

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

### Greg von Winckel (view profile)

03 Dec 2005 (Updated )

Compute Gauss nodes and weights for a triangle

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)
14 Jun 2012 Cássio

### Cássio (view profile)

Hey Greg, great code!
Thanks

01 Apr 2012 sfilop

### sfilop (view profile)

28 Jan 2011 Grzegorz Knor

### Grzegorz Knor (view profile)

Ok, it was my mistake, should be:
f=@(x,y) x.^2.*exp(y);
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);
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);
end

25 Nov 2010 Grzegorz Knor

### Grzegorz Knor (view profile)

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
?

Comment only
10 Feb 2010 saurabh gupta

### saurabh gupta (view profile)

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.

Comment only
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?

Comment only
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);

Comment only
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.

Comment only
06 Dec 2005 Greg von Winckel

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

Comment only
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.