File Exchange

image thumbnail

Gaussian Quadrature for Triangles

version 1.0 (1.84 KB) by

Compute Gauss nodes and weights for a triangle

9 Downloads

Updated

View License

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.

Comments and Ratings (16)

Cássio

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

sfilop

sfilop (view profile)

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

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
?

saurabh gupta

DAR .

Seems good. Too much simple thought.

meysam mousavi

Im Serious

Does not integrate correctly

YALA Hakim

I need this file for my code of FEM

zhenzhen Jin

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.

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?

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

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.

Greg von Winckel

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

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

Missing transpose symbol in comments.

MATLAB Release
MATLAB 7 (R14)

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

» Watch video