MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

### Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

# Thread Subject: Area integral over a triangle in MATLAB -- is numerical integration possible?

 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Vivek Saxena Date: 20 Mar, 2010 04:50:06 Message: 1 of 20 Hi, I'm wondering if it is possible to perform a double integral in MATLAB like \int \int f(x,y) dx dy where the domain of integration is a triangle, whose vertex coordinates are known? A full question is posted here: http://www.physicsforums.com/showthread.php?t=388101. Thanks in advance!
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Steven Lord Date: 20 Mar, 2010 05:01:13 Message: 2 of 20 "Vivek Saxena" wrote in message news:ho1k5u$hsc$1@fred.mathworks.com... > Hi, > > I'm wondering if it is possible to perform a double integral in MATLAB > like > > \int \int f(x,y) dx dy > > where the domain of integration is a triangle, whose vertex coordinates > are known? > > A full question is posted here: > http://www.physicsforums.com/showthread.php?t=388101. Look at QUAD2D. -- Steve Lord slord@mathworks.com comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Vivek Saxena Date: 20 Mar, 2010 05:35:05 Message: 3 of 20 "Steven Lord" wrote in message ... > > "Vivek Saxena" wrote in message > news:ho1k5u$hsc$1@fred.mathworks.com... > > Hi, > > > > I'm wondering if it is possible to perform a double integral in MATLAB > > like > > > > \int \int f(x,y) dx dy > > > > where the domain of integration is a triangle, whose vertex coordinates > > are known? > > > > A full question is posted here: > > http://www.physicsforums.com/showthread.php?t=388101. > > Look at QUAD2D. > > -- > Steve Lord > slord@mathworks.com > comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ > Thanks Steve. The approach I thought to define c(x) and d(x) is based on the equations of the sides of the triangle (defining the boundaries of the area element). But there are several possibilities depending on how the triangle is oriented with respect to the axes: an edge could have a zero slope, a nonzero finite slope or an infinite slope, and the different permutations of the ordering of the vertices will probably also play a role in defining c(x) and d(x) differently, for different triangles. Is there a neater way to employ QUAD2D to find the integral of f(x,y) over a triangular region with arbitrary vertex locations (x1, y1), (x2, y2) and (x3, y3)?
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Roger Stafford Date: 20 Mar, 2010 06:02:06 Message: 4 of 20 "Vivek Saxena" wrote in message ... > I'm wondering if it is possible to perform a double integral in MATLAB like > \int \int f(x,y) dx dy > where the domain of integration is a triangle, whose vertex coordinates are known?   If it happens you don't have quad2d, you can use the older dblquad which integrates over rectangular regions by making a change of variable from x,y to u,v defined by:  x = (1-u)*x2 + u*((1-v)*x2 + v*x3)  y = (1-u)*y2 + u*((1-v)*y2 + v*y3) where (x1,y1), (x2,y2), and (x3,y3) are the triangular region's three vertices. Of course you must change dx dy in the double integral to du dv with the proper Jacobian factor. What was the triangular region in x,y space then becomes a unit square in u,v space - the variables u and v each range from 0 to 1.   Of course it is more convenient to use quad2d if you have it. Roger Stafford
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Vivek Saxena Date: 20 Mar, 2010 06:29:07 Message: 5 of 20 "Roger Stafford" wrote in message ... > "Vivek Saxena" wrote in message ... > > I'm wondering if it is possible to perform a double integral in MATLAB like > > \int \int f(x,y) dx dy > > where the domain of integration is a triangle, whose vertex coordinates are known? > > If it happens you don't have quad2d, you can use the older dblquad which integrates over rectangular regions by making a change of variable from x,y to u,v defined by: > > x = (1-u)*x2 + u*((1-v)*x2 + v*x3) > y = (1-u)*y2 + u*((1-v)*y2 + v*y3) > > where (x1,y1), (x2,y2), and (x3,y3) are the triangular region's three vertices. Of course you must change dx dy in the double integral to du dv with the proper Jacobian factor. What was the triangular region in x,y space then becomes a unit square in u,v space - the variables u and v each range from 0 to 1. > > Of course it is more convenient to use quad2d if you have it. > > Roger Stafford Thanks Roger. But didn't you mean x = (1-u)*x1 + u*[(1-v)*x2 + v*x3] y = (1-u)*y1 + u*[(1-v)*y2 + v*y3] because the coordinate map you gave doesn't involve x1 and y1. But even so, this becomes a triangle in (u, v) space because u = 0, v = 0 => x = x1, y = y1 u = 0, v = 1 => x = x1, y = y1 u = 1, v = 0 => x = x2, y = y2 u = 1, v = 1 => x = x3, y = y3 Only the last 3 coordinates are meaningful essentially, because if u = 0, the way this map is constructed, the value of v becomes meaningless. Perhaps you meant something else?
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Vivek Saxena Date: 20 Mar, 2010 07:00:24 Message: 6 of 20 But I think this will do, because this transforms the triangle into a right angled triangle in (u, v) space with the vertices (0, 0), (0, 1), and (1,1). Now I just have to write the integrand as f(x(u,v), y(u,v)), multiply it by the Jacobian factor, and integrate from u = 0 to 1 and v = 0 to u (and not v = 0 to 1). Is it possible to do this symbolically? That is, what we now have to do is \int_{u=0}^{u=1} \int_{v=0}^{v=u} g(u, v) du dv What limits does one specify to quad2d here?
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Roger Stafford Date: 20 Mar, 2010 07:33:03 Message: 7 of 20 "Vivek Saxena" wrote in message ... > Thanks Roger. But didn't you mean > > x = (1-u)*x1 + u*[(1-v)*x2 + v*x3] > y = (1-u)*y1 + u*[(1-v)*y2 + v*y3] > > because the coordinate map you gave doesn't involve x1 and y1. But even so, this becomes a triangle in (u, v) space because > > u = 0, v = 0 => x = x1, y = y1 > u = 0, v = 1 => x = x1, y = y1 > u = 1, v = 0 => x = x2, y = y2 > u = 1, v = 1 => x = x3, y = y3 > > Only the last 3 coordinates are meaningful essentially, because if u = 0, the way this map is constructed, the value of v becomes meaningless. Perhaps you meant something else?   You are quite right in your first assertion, Vivek. It was a typing error and should be as you say, x = (1-u)*x1 + u*((1-v)*x2 + v*x3) y = (1-u)*y1 + u*((1-v)*y2 + v*y3)   However I disagree with your second statement. The interior of the triangle in x and y becomes the full interior of the unit square in u & v. At the far end the point (x1,y1) becomes stretched out to the line segment connecting (0,0) and (0,1) in u,v space. Just shy of (x1,y1) a very short strip in x,y space is stretched to the full unit length strip in u,v space. The transformation is valid because at the same time the Jacobian shrinks toward zero in compensation as you approach the point (x1,y1). Of course the inverse transformation from u,v to x,y becomes singular at that point but it doesn't matter in finding the double integral. It only means that dblquad will be a little wasteful in the number of points it uses toward the u = 0 end.   If you doubt the above statement, try this: x1 = randn; y1 = randn; % Random vertices for triangle in x,y x2 = randn; y2 = randn; x3 = randn; y3 = randn; N = 2000; u = rand(1,N); % u & v fill the unit square in u,v space v = rand(1,N); x = (1-u)*x1 + u.*((1-v)*x2 + v*x3); % x and y will fill the triangle y = (1-u)*y1 + u.*((1-v)*y2 + v*y3); plot([x1,x2,x3,x1],[y1,y2,y3,y1],'r-',x,y,'y.') axis equal You will see that the red triangle becomes nicely filled up with yellow dots even though they become tightly packed toward the (x1,y1) vertex, (where the Jacobian is getting smaller.) Roger Stafford
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Bruno Luong Date: 20 Mar, 2010 07:36:05 Message: 8 of 20 "Vivek Saxena" wrote in message ... > Hi, > > I'm wondering if it is possible to perform a double integral in MATLAB like > > \int \int f(x,y) dx dy > > where the domain of integration is a triangle, whose vertex coordinates are known? > > A full question is posted here: http://www.physicsforums.com/showthread.php?t=388101. It looks like a wave equation. Usually in FEM integration perform using Gaussian quadrature formula. On triangle, it involve barycentric coordinates. There are several schemes more or less complex and accurate. Very accurate integration such as QUAD2 is not necessary in FEM, but take a lot of time. Bruno
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Roger Stafford Date: 20 Mar, 2010 07:42:03 Message: 9 of 20 "Vivek Saxena" wrote in message ... > But I think this will do, because this transforms the triangle into a right angled triangle in (u, v) space with the vertices (0, 0), (0, 1), and (1,1). > > Now I just have to write the integrand as f(x(u,v), y(u,v)), multiply it by the Jacobian factor, and integrate from u = 0 to 1 and v = 0 to u (and not v = 0 to 1). Is it possible to do this symbolically? That is, what we now have to do is > > \int_{u=0}^{u=1} \int_{v=0}^{v=u} g(u, v) du dv > > What limits does one specify to quad2d here?   No Vivek, you should integrate over the full square u = 0 to u = 1 and v = 0 to v = 1. dblquad will not let you do otherwise. It is the full u,v square that the x,y triangle is mapped into by the given transformation. It does not matter that it has a singularity at the (x1,y1) vertex. Roger Stafford
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Vivek Saxena Date: 20 Mar, 2010 09:32:04 Message: 10 of 20 "Roger Stafford" wrote in message ... > No Vivek, you should integrate over the full square u = 0 to u = 1 and v = 0 to v = 1. dblquad will not let you do otherwise. It is the full u,v square that the x,y triangle is mapped into by the given transformation. It does not matter that it has a singularity at the (x1,y1) vertex. > > Roger Stafford Hi Roger, I think I see your point. But what if I use quad2d? So, if I understand you correctly, the triangle is mapped to the square {(u,v): 0<= u <= 1, 0 <= v <= 1} in (u,v) space.
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Roger Stafford Date: 20 Mar, 2010 18:10:05 Message: 11 of 20 "Vivek Saxena" wrote in message ... > Hi Roger, I think I see your point. But what if I use quad2d?   If you use quad2d, then you could use the original x,y coordinates directly without any transformation, though that is a bit awkward because in general you have to break the integration into two separate ranges which depend on how the three points of the triangle are oriented. The triangle's edges may not be aligned with either of the coordinate axes.   It might be more convenient to use an affine transformation that transforms your triangle into another triangle that is oriented along the coordinate directions so that separate integration ranges would not be required. For example, you could say  x = x1 + u*(x2-x1) + v*(x3-x1)  y = y1 + u*(y2-y1) + v*(y3-y1) This transforms the triangular region u >= 0, v >= 0, u+v <= 1 into your given x,y triangle. Thus your quad2d in the u,v domain could have simple limits: u from 0 to 1 and v from 0 to 1-u, and that is within quad2d's capabilities.   Of course you must still use the correct Jacobian factor in the integral if you do this transformation, but it is now just a constant quantity, namely the ratio between the two triangles' respective areas. > So, if I understand you correctly, the triangle is mapped to the square {(u,v): 0<= u <= 1, 0 <= v <= 1} in (u,v) space. Yes, that is correct (for the earlier transformation using dblquad.) Roger Stafford
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Vivek Saxena Date: 21 Mar, 2010 03:02:02 Message: 12 of 20 > > Yes, that is correct (for the earlier transformation using dblquad.) > > Roger Stafford Hi again Roger. I thought about this a bit, and I'm unable to convince myself that the nonlinear transformation you wrote really does map a triangle into a square. I guess you're right, but when I did the integral, the values all came out wrong (I have an analytical solution to the final problem with which I can compare). And since this is the only additional term in my calculation, I have no reason to suspect the rest of the code which has been working so far. Anyhow, I ran the following code clear, close all; clc; x1 = 0, x2 = 1, x3 = 1, y1 = 0, y2 = 0, y3 = 1; c = 1; for u=0:0.01:1     for v=0:0.01:1         x(c) = (1-u)*x1 + u*((1-v)*x2+v*x3);         y(c) = (1-u)*y1 + u*((1-v)*y2+v*y3);         c = c + 1;     end end scatter(x,y); That is I begin with a triangle which is right angled, with its vertices in (x, y) space as (0,0), (0,1) and (1,1) for simplicity. Now I use this transformation to construct x and y points for u and v ranging from 0 to 1. This gives me the coordinates of all the points inside the triangle (as required). But here, I have assumed that the corresponding region in (u, v) space is a square (because I allow u and v to run freely). So the forward transformation is a map from a square to a triangle. I am convinced about that. But how do I prove that the inverse transformation maps the triangle to the square {(u,v) : 0<=u<=1, 0<=v<=1}? Perhaps this is a trivial question, but I can't seem to get my head around it. I am using dblquad, and my functions look like this: function f = abcintegrand(u,v,a,b,c,k_0,x1,x2,x3,y1,y2,y3) J = sqrt(-1); argexp = (1-u).*x1 + u.*((1-v).*x2 + v.*x3); jacob = (-x1+(1-v).*x2+v.*x3).*(-u.*y2+u.*y3)-(-y1+(1-v).*y2+v.*y3).*(-u.*x2+u.*x3); p = (a + b*((1-u).*x1 + u.*((1-v).*x2 + v.*x3)) + c*((1-u).*y1 + u.*((1-v).*y2 + v.*y3))).*exp(-J*k_0*argexp); f = p.*jacob; FUNCTION CALL: fun = @(u,v) abcintegrand(u,v,a(i),b(i),c(i),k_0,x1,x2,x3,y1,y2,y3); b_integral(e,i) = (E_0*k_0*k_0*(eps_r(e)-1/mu_r(e))/(2*D(e)))*dblquad(fun,0,1,0,1);
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Vivek Saxena Date: 22 Mar, 2010 01:14:04 Message: 14 of 20 Hi Roger, thanks for the detailed reply. I am familiar with the theory of the method (being an engineering student myself :-) ). That is not the problem. I was wondering however, that given the forward transformation (x, y) --> (u, v) which clearly maps the triangle into the square (simply because u and v are now independent), can we show that (u,v) --> (x,y) is unique. Anyway since my last post, I managed to convince myself about this. That was just my way of trying to motivate the choice of the transformation. Also, when I said that the values of the integral being computed are turning out to be wrong, I was referring not to the transformation but rather to my usage of dblquad. I'm inclined to believe that something is wrong with my usage of that function. I should've been more explicit. (Thanks for your help btw.)
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Bruno Luong Date: 22 Mar, 2010 07:20:05 Message: 15 of 20 "Vivek Saxena" wrote in message ... > Hi Roger, thanks for the detailed reply. I am familiar with the theory of the method (being an engineering student myself :-) ). That is not the problem. I was wondering however, that given the forward transformation (x, y) --> (u, v) which clearly maps the triangle into the square (simply because u and v are now independent), can we show that (u,v) --> (x,y) is unique. Anyway since my last post, I managed to convince myself about this. That was just my way of trying to motivate the choice of the transformation. > > Also, when I said that the values of the integral being computed are turning out to be wrong, I was referring not to the transformation but rather to my usage of dblquad. I'm inclined to believe that something is wrong with my usage of that function. I should've been more explicit. (Thanks for your help btw.) One way to check the correctness of the parametrization is to integrate a polynomial P(x) on the triangle. If P(x) the polynomial of - order 0: the integral result is P(x1)*|T| - order 1; the integral result is (P(x1)+P(x2)+P(x3)) * |T|/3 etc... x1, x2, x3 are three corners, and |T| is the area of triangle Once the implementation gives correct result, you can then replace P(x) by your function. Bruno
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Roger Stafford Date: 22 Mar, 2010 15:19:03 Message: 16 of 20 "Vivek Saxena" wrote in message ... > ......... > Also, when I said that the values of the integral being computed are turning out to be wrong, I was referring not to the transformation but rather to my usage of dblquad. I'm inclined to believe that something is wrong with my usage of that function. ....... ------------   One thing I forgot to ask you, Vivek. Are your values of k_0 and the ranges of x in the triangle such that an excessive number of cycles occur in the trigonometric function, exp(-i*k_0*x)? If so, numerical integration is a correspondingly difficult process and dblquad, or possibly even quad2d, may very well end up confused and erroneous unless you place a very tight tolerance setting in them. This is true whether or not you have performed a transformation of the variables, and it is also true even for single integration.   Quadrature routines can be temperamental beasts when faced with rapidly fluctuating integrands. You would do well to read the description of quad2d at  http://www.mathworks.com/access/helpdesk/help/techdoc/ref/quad2d.html to get an understanding of these difficulties even if you are only using dblquad. Roger Stafford
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Bruno Luong Date: 22 Mar, 2010 15:31:04 Message: 17 of 20 "Roger Stafford" wrote in message ... > One thing I forgot to ask you, Vivek. Are your values of k_0 and the ranges of x in the triangle such that an excessive number of cycles occur in the trigonometric function, exp(-i*k_0*x)? Roger, I can answer this question. Usually not: when the FEM is correctly setup then the mesh size must be small compare to wave number, thus exp(-i*k_0*x) would not have many cycles on the triangle element. If it's not the case then there is no use to solve FEM because the discretized solution will have huge error. Bruno
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Vivek Saxena Date: 22 Mar, 2010 16:01:06 Message: 18 of 20 "Bruno Luong" wrote in message ... > "Roger Stafford" wrote in message ... > > > One thing I forgot to ask you, Vivek. Are your values of k_0 and the ranges of x in the triangle such that an excessive number of cycles occur in the trigonometric function, exp(-i*k_0*x)? > > Roger, I can answer this question. Usually not: when the FEM is correctly setup then the mesh size must be small compare to wave number, thus exp(-i*k_0*x) would not have many cycles on the triangle element. If it's not the case then there is no use to solve FEM because the discretized solution will have huge error. > > Bruno Hi Roger and Bruno. Thanks for your replies and help so far. I'm using a unit wavelength, and k_0 = 2*pi as a result. Also, the mesh size has been kept small enough -- actually this is a textbook problem but I have to solve it with a more complicated boundary condition than is given in the standard books. So these particular double integrals arise in the second order absorbing boundary condition, but not in the formulation with the first order absorbing boundary condition (for which the code runs perfectly). In fact, theory tells me that the second order ABC should work fine with a coarser discretization as compared to the first order ABC (if things are going fine). But I am being conservative and keeping the mesh size just as fine as it was earlier. Roger, I think it is likely that the program is plagued by some numerical errors, as you have pointed out. (I have gone through the code and ensured that there are no logical or syntactic errors.) I must admit that I do not know these techniques very well theoretically, but I will try and play with the options. If you have more suggestions on how the tolerances could be set in view of the kind of function that the integrand is, please do advise me.
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Roger Stafford Date: 23 Mar, 2010 17:03:23 Message: 19 of 20 "Bruno Luong" wrote in message ... > One way to check the correctness of the parametrization is to integrate a polynomial P(x) on the triangle. If P(x) the polynomial of > - order 0: the integral result is P(x1)*|T| > - order 1; the integral result is (P(x1)+P(x2)+P(x3)) * |T|/3 > etc... > > x1, x2, x3 are three corners, and |T| is the area of triangle > > Once the implementation gives correct result, you can then replace P(x) by your function. > > Bruno ------------   Hi Bruno. Your "etc..." intrigued me, so I worked out what the "order 2" test would be. I wasn't aware of this identity before. What would an order 3 test be?   The double integral of any quadratic polynomial,  P(x,y) = A*x^2+B*x*y+C*y^2+D*x+E*y+F , over a triangle is equal to the triangle's area, multiplied by one quarter of the polynomial's average value at the three vertices plus three quarters of its value at the triangle's centroid (located at the averages of the three respective coordinates.) Roger Stafford
 Subject: Area integral over a triangle in MATLAB -- is numerical integration possible? From: Bruno Luong Date: 23 Mar, 2010 18:11:05 Message: 20 of 20 "Roger Stafford" wrote in message ... > "Bruno Luong" wrote in message ... > > One way to check the correctness of the parametrization is to integrate a polynomial P(x) on the triangle. If P(x) the polynomial of > > - order 0: the integral result is P(x1)*|T| > > - order 1; the integral result is (P(x1)+P(x2)+P(x3)) * |T|/3 > > etc... > > > > x1, x2, x3 are three corners, and |T| is the area of triangle > > > > Once the implementation gives correct result, you can then replace P(x) by your function. > > > > Bruno > ------------ > Hi Bruno. Your "etc..." intrigued me, so I worked out what the "order 2" test would be. I wasn't aware of this identity before. What would an order 3 test be? > > The double integral of any quadratic polynomial, > > P(x,y) = A*x^2+B*x*y+C*y^2+D*x+E*y+F , > > over a triangle is equal to the triangle's area, multiplied by one quarter of the polynomial's average value at the three vertices plus three quarters of its value at the triangle's centroid (located at the averages of the three respective coordinates.) > > Roger Stafford Hi Roger, The order 3 polynomial, the formula is I = 9/20*T*[P((x1+x2+x3)/3)] + 2/15*T*[P((x1+x2)/2)+P((x2+x3)/2)+P((x3+x1)/2)] + 1/20*T*[P(x1) + P(x2) + P(x3)]. Those kind of Gauss's formula are derived mostly for FEM (I copy it from a book). Bruno

## Tags for this Thread

### What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.