Integrate a function over a triangle area

Hi, I'm trying to integrate a function over the area of a triangle, that is given by a set of 3 arbitrary points. My problem is that because the points are arbitrary, my program doesn't take into account that dA doesn't really equal dy.dx. Now, is there a predefined function that can do the integration over the area given by 3 points?

Here's my current code:

function [Ke]=Ke(F,ex,ey)
%F=@(x,y)(x/2 - y/2 + 1/2)^2 + (x/2 + y/2 - 1/2)^2 + y^2
%ex=[0 -1 1]
%ey=[1 0 0]

syms x y
x1= ex(1);
x2= ex(2);
x3= ex(3);
y1= ey(1);
y2= ey(2);
y3= ey(3);
Ke=zeros(3,3)
for i=1:3;
    for j=1:3;
       
        Ke(i,j)=integral2(F,x,ex(i),ex(j),ey(i),ey(j))
        K=K+Ke(i,j)
    end
end

[Bruno Luong]

IMO you have a much bigger problem than not knowing dA.

Your code adds the integrals of function 9 rectanges (3 x 3 pairs of x and y intervals), which won't give the integral on triangle at all.

This code is emply incorrect to start with.

[John D'Errico]

Is your function truly known, and as simple as the one you show? If so, you can write an analytical integral, there is no need for a call to integral2.

[Jafar Alrashdan]

The function isn't always the same as this one, but it will always be similar. So I need a general code that can integrate any function over any arbitrary set of points

Answer 1

[John D'Errico]

I'm not sure what the answer is, because there are multiple question left unresolved at this point. There are several approaches, thus, given a completely general function, and a list of points, you need first to answer what is the domain of integration. Can you evaluate the function at intermediate points? Is your function a simple one, as in the question? Or is that just some trivial example? Is your function well behaved, or can it get potentially nasty.

So, first, I suggest that you probably need to triangulate the domain. A list of points is not a triangulation. And this begs the question of whether your domain is a convex one, because tools like delaunay triangulations assume the domain is convex, because the resulting triangulation will always fill the convex hull of your point set.

Next, there are many ways to form an integration over a 2-d triangulated domain. For example, if you are unwilling to evaluate the function at new points inside each triangle, then at best you can use simple rules that presume the function is linear over each triangle. Thus a Guassian integration over a triangulated domain, where you have only function values available at the vertices of each triangle is REALLY easy to compute.

Moving upwards, if you can evaluate the function at new points inside each triangle, there are Gaussian integrations schemes over a triangle. The simplest ones involve an implicit transformation (a mapping) of the triangle into a square. Then use a tensor product mesh inside the square. With the correct transformations, this problem too is solvable. In fact, I have provided code for this. It may even be posted on the FEX. I can't recall.

Next, there are ways to write an adaptive scheme for integration on a triangulated mesh, where you use a sequence of nested integral estimates. The trick is to carefully re-use information, else this becomes inefficient, but to then apply Richardson extrapolation techniques to the nested sequence of integral estimates to gain a higher order estimate. I've written code for all this too, but adaptive 2-d integration over a general domain can be tricky to know when to stop refining the estimates based on a convergence tolerance. So, while I've written code for these things, it can be quite tricky to get just right.

Anyway, I'm not at all sure this is your real queation. Are you simply asking how to form a low order integral of a function over a general triangle, because you don't know what you need to do when the area of the triangle is not 1/2?

Or are you asking for something more general, thus how do you perform the integration over some general domain?

For example, you could simply use Greg von Winckel's nice tool, simplexquad. It is on the file exchange. It only applies to a simple simplex, however, we can do this:

xy = [ex;ey]'
xy =
     0     1
    -1     0
     1     0
F=@(x,y) (x/2 - y/2 + 1/2).^2 + (x/2 + y/2 - 1/2).^2 + y.^2;
[V,W]=simplexquad(2,xy)
V =
         -0.56991786405571         0.280019915499074
        0.0235077025419343         0.666390246014701
        -0.364929058779244        0.0750311102226081
         0.511339220293019         0.178558728263616
W =
         0.181958618256023
         0.318041381743977
         0.181958618256023
         0.318041381743977
W'*F(V(:,1),V(:,2))
ans =
                       0.5

And we see the integral of F over that domain, exact for the low order polynomial you have supplied, is 1/2. If you increase the value of n, it will yield the same result, as long as n>= 2 in the call to simplexquad.

Given a more complicated function over that domain however, we would see that you need to use a higher order Gaussian integration.

F = @(x,y) sin(x+y).*cos(x-2*y);
for n = 1:9
  [V,W]=simplexquad(n,xy);
  I = W'*F(V(:,1),V(:,2))
end
I =
         0.257138144005872
I =
          0.18489244716809
I =
         0.188773506173431
I =
          0.18871206014412
I =
         0.188712578059871
I =
          0.18871257529596
I =
         0.188712575306135
I =
         0.188712575306107
I =
         0.188712575306107

So it appears as if the Gaussian integration over the triangle has converged to double precision accuracy by the time we got to n==8, requiring 64 function evaluations. The problem is, I did not know that it has converged until I tried n=9 to see the result was the same.

The adaptive code I wrote is not as good, because it required more function evaluations. That is to be expected however.

options.plot = 1;
options.tol = 1.e-8;
I = simplxint(F,[1 2 3],,options)
I =
          0.18871257530507

The final triangulation arrived at was apparently:

to meet the integration tolerance as requested of 1e-8.

So again, I don't know exactly what it is that you are looking to do, or how I might help you. Simplest would just be download Greg'scode from the File Exchange. Eveneasier, I've attached it to this answer.

[Jafar Alrashdan]

simplexquad would solve this problem! thanks for the tip!

[John D'Errico]

Just in case you want the link, even though I did attach simplexquad, here are some useful links:

https://www.mathworks.com/matlabcentral/fileexchange/9435-n-dimensional-simplex-quadrature

https://www.mathworks.com/matlabcentral/fileexchange/9230-gaussian-quadrature-for-triangles

The first is simplexquad, the second link is one he posted for the specific case of a triangle.

Answer 2

[Bruno Luong]

Here is a method.

Assuming you triangle T has 3 vertexes P1, P2, P3, each vertex Pi has corrdinates (xi,yi).

Any point P inside a triangle can be written as barycentric coordinates (w1,w2,w3)

P = w1*P1 + w2*P2 + w3*P3

with wi such that

0 <= wi <=0

and

w1+w2+w3 = 1.

Integral of f on T written in barycentric cooordinates becomes

|T|/2 Integral f(w1*P1+w2*P2+(1-w1-w2)*P3) (dw1*dw2)

the integral is carried out on the "right" triangle domain T12 (where |T| is the area of the original triangle T given later):

T12: = { (w1,w2): 0<=w1<=1; 0<=w2<=1-w1 }.

So what you should do is compute this double-cascaded integrals:

|T|/2 * integral_(w1 in (0,1)) integral (w2 in (0,1-w1)) f(w1*P1+w2*P2+(1-w1-w2)*P3) dw2 dw1.

The area |T| (dA) can be computed

 |T| = 0.5* abs(det (P2-P1,P3-P1))

MATLAB code

function I = TriIntegral(f, Tx, Ty)
% I = TriIntegral(f, Tx, Ty)
% 2D integration of f on a triangle
% INPUTS:
%   - f is the vectorized function handle that when calling f(x,y) returns
%       function value at (x,y), x and y are column vectors
%   - Tx,Ty are is two vectors of length 3, coordinates of the triangle
% OUTPUT
%   I: integral of f in T

T = [Tx(:), Ty(:)];
I = integral2(@(s,t) fw(f,s,t,T),0,1,0,1);
A = det(T(2:3,:)-T(1,:));
I = I*abs(A);
end % TriIntegral

%%
function y = fw(f, s, t, T)
sz = size(s);
w1 = (1-s); % Bug fix
w2 = s.*t;
w3 = 1-w1-w2;
P = [w1(:),w2(:),w3(:)] * T;
y = feval(f,P(:,1),P(:,2));
y = s(:).*y(:);
y = reshape(y,sz);
end

Apply to your example:

F=@(x,y)(x/2 - y/2 + 1/2).^2 + (x/2 + y/2 - 1/2).^2 + y.^2;
ex=[0 -1 1;
ey=[1 0 0];
TriIntegral(F,ex,ey)

>> ans =

    0.5000

[Bruno Luong]

Note:

I = integral_(w1 in (0,1)) integral (w2 in (0,w1)) g(w1,w2) dw1 dw2

Can be transformed to an integral on square, thus using integral_2

I = integral_(s in (0,1)) integral (t in (0,1)) g(s,s*t)*s ds dt

[Bruno Luong]

Note 2:

If you have f polynomial function; the integration can be computed by using Gauss and weight points. Checkout https://en.wikipedia.org/wiki/Gaussian_quadrature

[Jafar Alrashdan]

Thank you for your detailed answer, I'll have to read about barycentric cooordinates.

As for Gauss, I've already explored that option, but found it's implementation into a trianglur region too complicated.

[Bruno Luong]

Not complicated at all if you don't want to reinvent the wheel.

You have many File Exchange available that compute the Points and Weights with respect to the Barycentric coordinates (again), and apply it is simply a matter of summing up the function.

[Jafar Alrashdan]

Thanks again, I'll look it up

[John D'Errico]

As I showed in my response, a Gaussian integration is simple, if you use Greg's code to give a Gaussian integration over a simplex.

[Bruno Luong]

wops, I made a mistake, w2 should be in (0,1-w1). Code and text corrected accordingly

[Bruno Luong]

Try the same test than John's, it differs by 2 last digits.

The advantage of using integral2 over Gaussian quadrature is that it will automatically adapt to the function smoothness and you don't have to ask whereas it converges or not.

F = @(x,y) sin(x+y).*cos(x-2*y)
ex=[0 -1 1];
ey=[1 0 0];
TriIntegral(F,ex,ey)

ans =

   0.188712575302678

Answer 3

[madhan ravi]

1) Find the equation of the 3 lines connecting the three points.

2)Split the domain according to upper limit and lower limit by substituting the points in the domain you can determine it.

3) Use double integration in each domain and sum up all the integrals to attain the final triangulare area.

Read Thomas calculus to know more

[Jafar Alrashdan]

I can easily do it analaytically if the function and points are the always known. But my goal is to just give an arbitrary function and set points, and the matlab function would integrate over it.

I'm not sure if my explanation is clear enough

[madhan ravi]

I know what you mean but the domain can only be defined by the user becuase matlab won't analayse the upper and lower limits by itself.

[Jafar Alrashdan]

Thanks for your help, I guess I have to figure out a way to program a function for it to know the limits

[madhan ravi]

Your welcome , if you do it publish it as file exchange

[Jafar Alrashdan]

I will