Gaussian Quadrature for Triangles

Compute Gauss nodes and weights for a triangle

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)

[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.

[Greg von Winckel]

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

[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.

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

[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?

[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.

[zhenzhen Jin]

 

[YALA Hakim]

I need this file for my code of FEM

[Im Serious]

Does not integrate correctly

[meysam mousavi]

 

[DAR .]

Seems good. Too much simple thought.

[saurabh gupta]

 

[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
?

[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