View License

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Highlights from
Gaussian Quadrature for Triangles

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

Gaussian Quadrature for Triangles



03 Dec 2005 (Updated )

Compute Gauss nodes and weights for a triangle

| Watch this File

File Information

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?

01 Apr 2012 sfilop

sfilop (view profile)

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

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]);

[X,Y,Wx,Wy]=triquad(10,[0 0; 0 2; 2 0]);

[X,Y,Wx,Wy]=triquad(20,[0 0; 0 2; 2 0]);

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

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

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:


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

07 Dec 2005

Missing transpose symbol in comments.

Contact us