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.