Discover MakerZone

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

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
How to integration a polynomial on a polygon?

Subject: How to integration a polynomial on a polygon?

From: zedong

Date: 21 Dec, 2008 04:28:03

Message: 1 of 16

How to integration a polynomial on a polygon?
for example:In 2 dimension:
I have a function f(x,y)=2*x+6*y+7
I want to get the integration over a polygon domain:
for example:
[1 1;0.4 0.4;0 1;0 0]
is a polygon boundary coordinates.
Remark:
(1) The polygon may not be convex.
(2) I will have to computation for many times of different function and different polygon.

Subject: How to integration a polynomial on a polygon?

From: zedong

Date: 21 Dec, 2008 04:33:02

Message: 2 of 16

Could you give me a general algorithm?Or is there any easy implementation?
Thank you very much.

Subject: How to integration a polynomial on a polygon?

From: Walter Roberson

Date: 21 Dec, 2008 05:06:50

Message: 3 of 16

zedong wrote:

> How to integration a polynomial on a polygon?
> for example:In 2 dimension:
> I have a function f(x,y)=2*x+6*y+7
> I want to get the integration over a polygon domain:
> for example:
> [1 1;0.4 0.4;0 1;0 0]
> is a polygon boundary coordinates.
> Remark:
> (1) The polygon may not be convex.
> (2) I will have to computation for many times of different function and different polygon.

> Could you give me a general algorithm?Or is there any easy implementation?

Rotate the coordinate space around one of the line segment endpoints so that the
other line segment endpoint was on the y axis. The problem of finding the
intersection of the polynomial with the line segment is then reduced to the
problem of finding the zeros of a polynomial of the same degree (since
rotation around a point is a linear transformation), and you already know
the general algorithm and easy implementation for that, right?

--
.signature note: I am now avoiding replying to unclear or ambiguous postings.
Please review questions before posting them. Be specific. Use examples of what you mean,
of what you don't mean. Specify boundary conditions, and data classes and value
relationships -- what if we scrambled your data or used -Inf, NaN, or complex(rand,rand)?

Subject: How to integration a polynomial on a polygon?

From: Roger Stafford

Date: 21 Dec, 2008 09:56:04

Message: 4 of 16

"zedong
" <zdongwu@gmail.com> wrote in message <gikgkj$4qp$1@fred.mathworks.com>...
> How to integration a polynomial on a polygon?
> for example:In 2 dimension:
> I have a function f(x,y)=2*x+6*y+7
> I want to get the integration over a polygon domain:
> for example:
> [1 1;0.4 0.4;0 1;0 0]
> is a polygon boundary coordinates.
> Remark:
> (1) The polygon may not be convex.
> (2) I will have to computation for many times of different function and different polygon.

  For linear functions such as in your f(x,y), the integral over a polygonal domain can be evaluated by a comparatively simple summation. I'll just give the method for f(x,y) = x.

  Let x and y be vectors of the polygon vertices' coordinates. Then do:

 x1 = x([2:end,1]);
 y1 = y([2:end,1]);
 I = sum((x.^2+x.*x1+x1.^2).*(y1-y))/6;

I is the desired integral of f(x,y) = x over the polygonal area.

  There are similar formulas for f(x,y) = y and f(x,y) = 1 (area), and you could combine them in a linear combination to evaluate the function above.

  Note that it is assumed that the polygon vertices in the vectors are sequenced in counterclockwise order around the polygon. Also it is assumed that the polygon boundary never crosses itself, as for example in a figure eight shape. The polygon does not have to be convex.

  Such formulas become considerably more complicated for functions like x^2, x*y, x^2*y, etc. or linear combinations of these, but they are possible.

Roger Stafford

Subject: How to integration a polynomial on a polygon?

From: John D'Errico

Date: 21 Dec, 2008 11:27:02

Message: 5 of 16

"zedong
" <zdongwu@gmail.com> wrote in message <gikgkj$4qp$1@fred.mathworks.com>...
> How to integration a polynomial on a polygon?
> for example:In 2 dimension:
> I have a function f(x,y)=2*x+6*y+7
> I want to get the integration over a polygon domain:
> for example:
> [1 1;0.4 0.4;0 1;0 0]
> is a polygon boundary coordinates.
> Remark:
> (1) The polygon may not be convex.
> (2) I will have to computation for many times of different function and different polygon.

You asked for a general algorithm, an "easy"
implementation. I can give you an algorithm
easily enough. Is it easy? That is probably a
function of your skill at matlab. While I have
written code to do this, I've not yet put it on
the file exchange.

1. Triangulate the polygon.

2. Integrate your function over each triangular
piece.

The first step will require one of several
techniques to triangulate a valid polygon.
You will find mesh2d on the file exchange.
So step 1 is done for you.

I've already put a function called tessquad
on the file exchange that will do a gaussian
integration of a function on a triangulated
domain. So step 2 is also done for you.

This scheme will work on ANY polygon as long
as it does not cross itself, convex or not.

HTH,
John

Subject: How to integration a polynomial on a polygon?

From: Roger Stafford

Date: 21 Dec, 2008 23:11:02

Message: 6 of 16

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gil3rj$ss6$1@fred.mathworks.com>...
> For linear functions such as in your f(x,y), the integral over a polygonal domain can be evaluated by a comparatively simple summation. I'll just give the method for f(x,y) = x.
>
> Let x and y be vectors of the polygon vertices' coordinates. Then do:
>
> x1 = x([2:end,1]);
> y1 = y([2:end,1]);
> I = sum((x.^2+x.*x1+x1.^2).*(y1-y))/6;
>
> I is the desired integral of f(x,y) = x over the polygonal area.
>
> There are similar formulas for f(x,y) = y and f(x,y) = 1 (area), and you could combine them in a linear combination to evaluate the function above.
>
> Note that it is assumed that the polygon vertices in the vectors are sequenced in counterclockwise order around the polygon. Also it is assumed that the polygon boundary never crosses itself, as for example in a figure eight shape. The polygon does not have to be convex.
>
> Such formulas become considerably more complicated for functions like x^2, x*y, x^2*y, etc. or linear combinations of these, but they are possible.
>
> Roger Stafford

  The following allows you to generalize that formula I gave you, Zedong. Unfortunately you will have to convert its format to that of up-to-date symbolic toolboxes. Mine is an older version. Suppose you want to integrate the function f(x,y) = x*y. Then do this:

 a = 1; b = 1;
 f = subs(subs('x^a*y^b',a,'a'),b,'b');
 I1 = int(int(f,'y',0,'y2/x2*x'),'x',0,'x2');
 I2 = int(int(f,'y',0,'y1/x1*x'),'x',0,'x1');
 I3 = int(int(f,'y',0,'(y1*x2-y2*x1+(y2-y1)*x)/(x2-x1)'),'x','x2','x1');
 I = factor(simplify(symop(I1,'-',I2,'+',I3)));
 pretty(I)
 
       1/24 (- y1 x2 + y2 x1) (2 x1 y1 + y2 x1 + y1 x2 + 2 x2 y2)

Corresponding expression for I:

 I = '(y2*x1-y1*x2)*(y1*(2*x1+x2)+y2*(2*x2+x1))/24';

  Now substitute x for x1, x1 for x2, y for y1, and y1 for y2 in this expression and use the method I described earlier with polygon vectors x and y.

 x1 = x([2:end,1]);
 y1 = y([2:end,1]);
 K = sum((y1.*x-y.*x1).*(y.*(2.*x+x1)+y1.*(2.*x1+x)))/24;

K is your integral.

  Note 1: "symop(I1,'-',I2,'+',I3)" is just my symbolic toolbox's crude way of combining the expressions: I1-I2+I3.

  Note 2: The above expression for I is the integral of f(x,y) = x*y over the area of the triangle (0,0), (x1,y1), (x2,y2). In the summation all areas of integration lying outside the polygon eventually cancel out.

  Note 3: At about degree 5, these formulas begin to be uncomfortably large and cumbersome so you might prefer numerical methods for polynomials of that high a degree.

Roger Stafford

Subject: How to integration a polynomial on a polygon?

From: zedong

Date: 22 Dec, 2008 01:01:53

Message: 7 of 16

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gimie6$gr7$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gil3rj$ss6$1@fred.mathworks.com>...
> > For linear functions such as in your f(x,y), the integral over a polygonal domain can be evaluated by a comparatively simple summation. I'll just give the method for f(x,y) = x.
> >
> > Let x and y be vectors of the polygon vertices' coordinates. Then do:
> >
> > x1 = x([2:end,1]);
> > y1 = y([2:end,1]);
> > I = sum((x.^2+x.*x1+x1.^2).*(y1-y))/6;
> >
> > I is the desired integral of f(x,y) = x over the polygonal area.
> >
> > There are similar formulas for f(x,y) = y and f(x,y) = 1 (area), and you could combine them in a linear combination to evaluate the function above.
> >
> > Note that it is assumed that the polygon vertices in the vectors are sequenced in counterclockwise order around the polygon. Also it is assumed that the polygon boundary never crosses itself, as for example in a figure eight shape. The polygon does not have to be convex.
> >
> > Such formulas become considerably more complicated for functions like x^2, x*y, x^2*y, etc. or linear combinations of these, but they are possible.
> >
> > Roger Stafford
>
> The following allows you to generalize that formula I gave you, Zedong. Unfortunately you will have to convert its format to that of up-to-date symbolic toolboxes. Mine is an older version. Suppose you want to integrate the function f(x,y) = x*y. Then do this:
>
> a = 1; b = 1;
> f = subs(subs('x^a*y^b',a,'a'),b,'b');
> I1 = int(int(f,'y',0,'y2/x2*x'),'x',0,'x2');
> I2 = int(int(f,'y',0,'y1/x1*x'),'x',0,'x1');
> I3 = int(int(f,'y',0,'(y1*x2-y2*x1+(y2-y1)*x)/(x2-x1)'),'x','x2','x1');
> I = factor(simplify(symop(I1,'-',I2,'+',I3)));
> pretty(I)
>
> 1/24 (- y1 x2 + y2 x1) (2 x1 y1 + y2 x1 + y1 x2 + 2 x2 y2)
>
> Corresponding expression for I:
>
> I = '(y2*x1-y1*x2)*(y1*(2*x1+x2)+y2*(2*x2+x1))/24';
>
> Now substitute x for x1, x1 for x2, y for y1, and y1 for y2 in this expression and use the method I described earlier with polygon vectors x and y.
>
> x1 = x([2:end,1]);
> y1 = y([2:end,1]);
> K = sum((y1.*x-y.*x1).*(y.*(2.*x+x1)+y1.*(2.*x1+x)))/24;
>
> K is your integral.
>
> Note 1: "symop(I1,'-',I2,'+',I3)" is just my symbolic toolbox's crude way of combining the expressions: I1-I2+I3.
>
> Note 2: The above expression for I is the integral of f(x,y) = x*y over the area of the triangle (0,0), (x1,y1), (x2,y2). In the summation all areas of integration lying outside the polygon eventually cancel out.
>
> Note 3: At about degree 5, these formulas begin to be uncomfortably large and cumbersome so you might prefer numerical methods for polynomials of that high a degree.
>
> Roger Stafford



Thank you very much for your kind help.I am learning your algorithm.Thank you very much.It's very beautiful.

Subject: How to integration a polynomial on a polygon?

From: zedong

Date: 22 Dec, 2008 01:20:05

Message: 8 of 16

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gimie6$gr7$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gil3rj$ss6$1@fred.mathworks.com>...
> > For linear functions such as in your f(x,y), the integral over a polygonal domain can be evaluated by a comparatively simple summation. I'll just give the method for f(x,y) = x.
> >
> > Let x and y be vectors of the polygon vertices' coordinates. Then do:
> >
> > x1 = x([2:end,1]);
> > y1 = y([2:end,1]);
> > I = sum((x.^2+x.*x1+x1.^2).*(y1-y))/6;
> >
> > I is the desired integral of f(x,y) = x over the polygonal area.
> >
> > There are similar formulas for f(x,y) = y and f(x,y) = 1 (area), and you could combine them in a linear combination to evaluate the function above.
> >
> > Note that it is assumed that the polygon vertices in the vectors are sequenced in counterclockwise order around the polygon. Also it is assumed that the polygon boundary never crosses itself, as for example in a figure eight shape. The polygon does not have to be convex.
> >
> > Such formulas become considerably more complicated for functions like x^2, x*y, x^2*y, etc. or linear combinations of these, but they are possible.
> >
> > Roger Stafford
>
> The following allows you to generalize that formula I gave you, Zedong. Unfortunately you will have to convert its format to that of up-to-date symbolic toolboxes. Mine is an older version. Suppose you want to integrate the function f(x,y) = x*y. Then do this:
>
> a = 1; b = 1;
> f = subs(subs('x^a*y^b',a,'a'),b,'b');
> I1 = int(int(f,'y',0,'y2/x2*x'),'x',0,'x2');
> I2 = int(int(f,'y',0,'y1/x1*x'),'x',0,'x1');
> I3 = int(int(f,'y',0,'(y1*x2-y2*x1+(y2-y1)*x)/(x2-x1)'),'x','x2','x1');
> I = factor(simplify(symop(I1,'-',I2,'+',I3)));
> pretty(I)
>
> 1/24 (- y1 x2 + y2 x1) (2 x1 y1 + y2 x1 + y1 x2 + 2 x2 y2)
>
> Corresponding expression for I:
>
> I = '(y2*x1-y1*x2)*(y1*(2*x1+x2)+y2*(2*x2+x1))/24';
>
> Now substitute x for x1, x1 for x2, y for y1, and y1 for y2 in this expression and use the method I described earlier with polygon vectors x and y.
>
> x1 = x([2:end,1]);
> y1 = y([2:end,1]);
> K = sum((y1.*x-y.*x1).*(y.*(2.*x+x1)+y1.*(2.*x1+x)))/24;
>
> K is your integral.
>
> Note 1: "symop(I1,'-',I2,'+',I3)" is just my symbolic toolbox's crude way of combining the expressions: I1-I2+I3.
>
> Note 2: The above expression for I is the integral of f(x,y) = x*y over the area of the triangle (0,0), (x1,y1), (x2,y2). In the summation all areas of integration lying outside the polygon eventually cancel out.
>
> Note 3: At about degree 5, these formulas begin to be uncomfortably large and cumbersome so you might prefer numerical methods for polynomials of that high a degree.
>
> Roger Stafford

I am sorry.I think your method must be very beautiful.But I can't follow your algorithm.Could you tell me why this kind of combination summation is right?Thank you very much.If it's hard to explain, Could you tell me where can I find the theory?Thank you.

Subject: How to integration a polynomial on a polygon?

From: Roger Stafford

Date: 22 Dec, 2008 03:20:05

Message: 9 of 16

"zedong" <zdongwu@gmail.com> wrote in message <gimq05$2gv$1@fred.mathworks.com>...
> .......
> I am sorry.I think your method must be very beautiful. But I can't follow your algorithm. Could you tell me why this kind of combination summation is right? Thank you very much. If it's hard to explain, Could you tell me where can I find the theory?Thank you.
----------
  In justifying this method there are three distinct ideas to explain. The first is why the expressions I1, I2, and I3 yield the integral of f(x,y) over the triangular region (0,0), (x1,y1), (x2,y2). The second is why the sum of all such triangular regions will yield the integration taken over only the area of the polygon. The third is why the summation code works.

A. The Triangle

  Suppose the triangle's vertices are O = (0,0), P = (x1,y1) = (3,1), and Q = (x2,y2) = (1,2). The line OQ has the equation y = y2/x2*x. If we project Q perpendicularly down to the x-axis at point R = (1,0) we have a right triangle ORQ. The quantity int(f,'y',0,'y2/x2*x') integrates f(x,y) with respect to y along a vertical line at x until meeting the line y = y2/x2*x, and then int(int(f,'y',0,'y2/x2*x'),'x',0,'x2') integrates that answer with respect to x from x = 0 to x = x2. Hence we have the integral of f(x,y) over triangle ORQ. In a similar way I3 integrates f(x,y) over the trapezoid RSPQ where S is P projected orthogonally to the x-axis. Finally I2 is the integral over the right triangle OSP. When we form I = I1-I2+I3 we have the integral over just the triangle OPQ. The area outside triangle OPQ has been subtracted out with the -I2 step.

B. The Summation

  Now suppose our polygon consists of three points A = (1,-1), B = (1,0), C = (2,2). Using I above on triangle OAB gives the integral over that triangle. Using it on triangle OBC gives the integral over that triangle. In both cases the ordering was counterclockwise around the triangles. When we finally take the integral over triangle OCA, we are going clockwise instead of counterclockwise, so the I in that case is the negative of the integral of f(x,y) over triangle OCA. This has the effect of subtracting the area that had been encompassed in the first two integrations that lay outside triangle ABC. Hence we end up with just the integral over ABC.

  Of course there are very many combinations of arrangement of the polygon to be considered but you will find that in every case the areas covered outside the polygon are eventually cancelled out when the summation has been completed. You will note that we could have chosen some other point rather than the coordinate origin to base our triangles on. Its choice is arbitrary and we would still end up with the same answer using other base points. However, using the origin tends to make the formulas more simple looking than other choices.

C. The Code

  In writing x1 = x([2:end,1]) and y1 = y([2:end,1]) this places the x1 and y1 vectors one step ahead of the points given by x and y. When we are using {x(i),y(i)} together with {x1(i),y1(i)), we are combining two points with the origin (0,0) to find the integral over that triangle. The first point plays the role of (x1,y1) and the second the role of (x2,y2) in the formula. The last point (x(n),y(n)) is combined with the first point (x1(n),y1(n)) = (x(1),y(1)) which is just what is needed to go completely around the polygon. The 'sum' operator carries out the desired summation process.


  The nice thing about using the Symbolic Toolbox for this purpose is that any desired polynomial term x^a*y^b can be used in the integration, though things get more and more complicated when a and b are large.

  As you may have noted, the results here are exact, since the integration results in 'int' above are exact. There are no numerical quadrature approximations involved, except of course for the inevitable round off errors from carrying out the arithmetical operations specified by the formulas.

  Incidentally note that the formula I gave you earlier for the case f(x,y) = x is different from the one here. That is because the earlier one was based on trapezoids off to the left of the various polygonal segments rather than triangles projected to the origin. The quantities added are different but their sum is still remains the same.

Roger Stafford

Subject: How to integration a polynomial on a polygon?

From: Peter

Date: 22 Dec, 2008 06:06:37

Message: 10 of 16

On Dec 20, 8:33=A0pm, "zedong " <zdon...@gmail.com> wrote:
> Could you give me a general algorithm?Or is there any easy implementation=
?
> Thank you very much.

Here is another way to approach the problem that allows a systematic
formulation of a computer program for an arbitrary simple polygon and
an arbitrary polynomial containing a finite number of coefficients...

Suppose the polynomial you want to integrate is p(x,y). (Example: p
(x,y) =3D 1 + x +2xy + y^2).
Compute the polynomials P(x,y) and Q(x,y) such that (d/dy) P(x,y) =3D
(-1/2) * p(x,y) and (d/dy) Q(x,y) =3D (1/2) * p(x,y). Here (d/dx) means
the partial derivative with respect to x, and similarly for (d/dy).
For our example polynomial these are P(x,y) =3D (-1/6)(3y + 3xy + 3xy^2
+ y^3) and Q(x,y) =3D (1/4)(2x + x^2 + 2yx^2 + 2xy^2). Then the surface
integral over the polygon of p(x,y) dx dy is equal to the surface
integral over the polygon of [(d/dx) Q(x,y) - (d/dy) P(x,y)] dx dy
which is equal (by Green's theorem in the plane) to the closed line
integral of [P(x,y) dx + Q(x,y) dy] around the perimeter of the
polygon in the counter-clockwise direction. This latter integral can
be written as the sum of the line integrals along each side of the
polygon. Since each of these is a portion of a straight line, it is
trivially easy to parameterize and evaluate these integrals (the
integrand being merely a polynomial, whose anti-derivative can be
written down immediately).

I've used a similar approach in the past to evaluate the Fourier
transform of a polygonal functions defined on polygonal domains. The
resulting formulas are compact and efficient.

Hope this helps,
--Peter

Subject: How to integration a polynomial on a polygon?

From: Peter

Date: 22 Dec, 2008 06:18:49

Message: 11 of 16

On Dec 21, 10:06=A0pm, Peter <petersamsim...@hotmail.com> wrote:
> On Dec 20, 8:33=A0pm, "zedong =A0" <zdon...@gmail.com> wrote:
>
> > Could you give me a general algorithm?Or is there any easy implementati=
on?
> > Thank you very much.
>
> Here is another way to approach the problem that allows a systematic
> formulation of a computer program for an arbitrary simple polygon and
> an arbitrary polynomial containing a finite number of coefficients...
>
> Suppose the polynomial you want to integrate is p(x,y). =A0(Example: p
> (x,y) =3D 1 + x +2xy + y^2).
> Compute the polynomials P(x,y) and Q(x,y) such that (d/dy) P(x,y) =3D
> (-1/2) * p(x,y) and (d/dy) Q(x,y) =3D (1/2) * p(x,y). Here (d/dx) means
> the partial derivative with respect to x, and similarly for (d/dy).
> For our example polynomial these are P(x,y) =3D (-1/6)(3y + 3xy + 3xy^2
> + y^3) and Q(x,y) =3D (1/4)(2x + x^2 + 2yx^2 + 2xy^2). =A0Then the surfac=
e
> integral over the polygon of p(x,y) dx dy is equal to the surface
> integral over the polygon of [(d/dx) Q(x,y) - (d/dy) P(x,y)] dx dy
> which is equal (by Green's theorem in the plane) to the closed line
> integral of [P(x,y) dx + Q(x,y) dy] around the perimeter of the
> polygon in the counter-clockwise direction. =A0This latter integral can
> be written as the sum of the line integrals along each side of the
> polygon. =A0Since each of these is a portion of a straight line, it is
> trivially easy to parameterize and evaluate these integrals (the
> integrand being merely a polynomial, whose anti-derivative can be
> written down immediately).
>
> I've used a similar approach in the past to evaluate the Fourier
> transform of a polygonal functions defined on polygonal domains. =A0The
> resulting formulas are compact and efficient.
>
> Hope this helps,
> --Peter

After a few more minutes of thought, I realized that you don't need
both P and Q. For instance, if you define Q(x,y) so that (d/dx) Q
(x,y) =3D p(x,y), you can simply choose P(x,y) =3D 0, and you only have a
single function to integrate around the perimeter of the polygon.

--Peter

Subject: How to integration a polynomial on a polygon?

From: Roger Stafford

Date: 22 Dec, 2008 08:14:06

Message: 12 of 16

Peter <petersamsimon2@hotmail.com> wrote in message <bb03a633-a402-444a-aed7-00ca729eba07@s16g2000vbp.googlegroups.com>...
> After a few more minutes of thought, I realized that you don't need
> both P and Q. For instance, if you define Q(x,y) so that (d/dx) Q
> (x,y) = p(x,y), you can simply choose P(x,y) = 0, and you only have a
> single function to integrate around the perimeter of the polygon.
>
> --Peter

  Using the d/dx Q(x,y) = p(x,y) and P(x,y) = 0 method would amount to evaluating the area integral by doing an iterated integration, first with respect to the x-direction and then with respect to the y-direction. The first integration results in the difference between the Q values at opposite ends of horizontal lines through the area. The second integration corresponds to the line integral of Q(x,y)dy along the left and right sides of the area.

  Yes, one could certainly do the problem this way. The first integration is very easily done for polynomials. However, I envision some difficulties with the complexity of the second integration. If, say, p(x,y) = x^3*y^2, with Q(x,y) = x^4*y^2/4, and we were then integrating between the points (x1,y1) and (x2,y2) as a line integral along the line segment, we would be dealing with an integral something like this

 int('y^2*((x1*y2-x2*y1+(x2-x1)*y)/(y2-y1))^4/4','y')

and I think this is where we would run into complexities of expression again. Of course, just I did above, we could chicken out and let the Symbolic Toolbox do all the hard work, but the resulting expression in terms of the various pairs of points (x1,y1) and (x2,y2) could still become rather annoyingly complex. It's not clear to me which method would lead to the greatest computing complexity. I would be interested in knowing the answer.

  Of course it must hold true that with either method after all integration is done. the integral values will necessarily both be the same function of the various points' coordinates, and I think that can be quite inherently messy for the higher degree polynomials. This seems to be the penalty one pays for obtaining exact answers as against numerical quadrature approximations.

Roger Stafford

Subject: How to integration a polynomial on a polygon?

From: zedong

Date: 22 Dec, 2008 08:14:06

Message: 13 of 16

Peter <petersamsimon2@hotmail.com> wrote in message <bb03a633-a402-444a-aed7-00ca729eba07@s16g2000vbp.googlegroups.com>...
> On Dec 21, 10:06=A0pm, Peter <petersamsim...@hotmail.com> wrote:
> > On Dec 20, 8:33=A0pm, "zedong =A0" <zdon...@gmail.com> wrote:
> >
> > > Could you give me a general algorithm?Or is there any easy implementati=
> on?
> > > Thank you very much.
> >
> > Here is another way to approach the problem that allows a systematic
> > formulation of a computer program for an arbitrary simple polygon and
> > an arbitrary polynomial containing a finite number of coefficients...
> >
> > Suppose the polynomial you want to integrate is p(x,y). =A0(Example: p
> > (x,y) =3D 1 + x +2xy + y^2).
> > Compute the polynomials P(x,y) and Q(x,y) such that (d/dy) P(x,y) =3D
> > (-1/2) * p(x,y) and (d/dy) Q(x,y) =3D (1/2) * p(x,y). Here (d/dx) means
> > the partial derivative with respect to x, and similarly for (d/dy).
> > For our example polynomial these are P(x,y) =3D (-1/6)(3y + 3xy + 3xy^2
> > + y^3) and Q(x,y) =3D (1/4)(2x + x^2 + 2yx^2 + 2xy^2). =A0Then the surfac=
> e
> > integral over the polygon of p(x,y) dx dy is equal to the surface
> > integral over the polygon of [(d/dx) Q(x,y) - (d/dy) P(x,y)] dx dy
> > which is equal (by Green's theorem in the plane) to the closed line
> > integral of [P(x,y) dx + Q(x,y) dy] around the perimeter of the
> > polygon in the counter-clockwise direction. =A0This latter integral can
> > be written as the sum of the line integrals along each side of the
> > polygon. =A0Since each of these is a portion of a straight line, it is
> > trivially easy to parameterize and evaluate these integrals (the
> > integrand being merely a polynomial, whose anti-derivative can be
> > written down immediately).
> >
> > I've used a similar approach in the past to evaluate the Fourier
> > transform of a polygonal functions defined on polygonal domains. =A0The
> > resulting formulas are compact and efficient.
> >
> > Hope this helps,
> > --Peter
>
> After a few more minutes of thought, I realized that you don't need
> both P and Q. For instance, if you define Q(x,y) so that (d/dx) Q
> (x,y) =3D p(x,y), you can simply choose P(x,y) =3D 0, and you only have a
> single function to integrate around the perimeter of the polygon.
>
> --Peter


Thanks you all.I must give my thanks to everybody for paying attention to this message.I can understand it now very clearly.Especially thanks should be given to Roger Stafford and Peter.Thank you for your kind help.

Subject: How to integration a polynomial on a polygon?

From: Roger Stafford

Date: 22 Dec, 2008 19:08:02

Message: 14 of 16

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message >
> .....
> int('y^2*((x1*y2-x2*y1+(x2-x1)*y)/(y2-y1))^4/4','y')
> .....

  Peter, when I wrote that expression for the line integral of Q(x,y) = x^4*y^2/4 from (x1,y1) to (x2,y2) I had in mind still obtaining a closed expression in terms of these points' coordinates for use in an exact summation. As a problem in numerical integration that would not be the right way to proceed. For example, it would blow up if y1 were equal to y2. It would be better to write something like this:

 int('Q(x,y)*sin(a12)','s',0,'s12')

for integration along each segment where s is the length along the segment with

 s12 = sqrt((x2-x1)^2+(y2-y1)^2)
 a12 = atan2(y2-y1,x2-x1)
 x = x1 + s*cos(a12)
 y = y1 + s*sin(a12)
 dy = sin(a12)*ds

  An analogous numerical line integration procedure exists for polynomials by computing their area integrals in terms of polygon segments projected triangularly to the origin. As in the case of your Q(x,y), the result of the first integration along the projecting lines to the origin is an easily-found explicit expression. The second integration along the polygonal line segments at the far end of these triangles could then be done numerically.

Roger Stafford

Subject: How to integration a polynomial on a polygon?

From: Peter

Date: 22 Dec, 2008 20:26:05

Message: 15 of 16

On Dec 22, 11:08=A0am, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> "Roger Stafford" <ellieandrogerxy...@mindspring.com.invalid> wrote in mes=
sage >
> > .....
> > =A0int('y^2*((x1*y2-x2*y1+(x2-x1)*y)/(y2-y1))^4/4','y')
> > .....
>
> =A0 Peter, when I wrote that expression for the line integral of Q(x,y) =
=3D x^4*y^2/4 from (x1,y1) to (x2,y2) I had in mind still obtaining a close=
d expression in terms of these points' coordinates for use in an exact summ=
ation. =A0As a problem in numerical integration that would not be the right=
 way to proceed. =A0For example, it would blow up if y1 were equal to y2. =
=A0It would be better to write something like this:
>
> =A0int('Q(x,y)*sin(a12)','s',0,'s12')
>
> for integration along each segment where s is the length along the segmen=
t with
>
> =A0s12 =3D sqrt((x2-x1)^2+(y2-y1)^2)
> =A0a12 =3D atan2(y2-y1,x2-x1)H
> =A0x =3D x1 + s*cos(a12)
> =A0y =3D y1 + s*sin(a12)
> =A0dy =3D sin(a12)*ds
>
> =A0 An analogous numerical line integration procedure exists for polynomi=
als by computing their area integrals in terms of polygon segments projecte=
d triangularly to the origin. =A0As in the case of your Q(x,y), the result =
of the first integration along the projecting lines to the origin is an eas=
ily-found explicit expression. =A0The second integration along the polygona=
l line segments at the far end of these triangles could then be done numeri=
cally.
>
> Roger Stafford

Hi, Roger.

I understood the original question to be asking for a numerical
evaluation of the integral, not a closed-form symbolic solution. The
Q polynomial coefficients could be calculated automatically by the
user's own code (no toolbox required), and their contributions along
each edge would be summed numerically inside of a loop over the edges,
again without requiring any external toolbox. Certainly things get a
lot more complicated if you want an explicit formula derived in terms
of the polygon vertex locations. I guess I was misleading when I
referred to the "formulas", by which I meant the formulas (finite
sums) needed for obtaining a numerical result.

--Peter

Subject: How to integration a polynomial on a polygon?

From: Roger Stafford

Date: 23 Dec, 2008 06:00:08

Message: 16 of 16

"zedong" <zdongwu@gmail.com> wrote in message <gini8e$3mo$1@fred.mathworks.com>...
> Thanks you all.I must give my thanks to everybody for paying attention to this message.I can understand it now very clearly.Especially thanks should be given to Roger Stafford and Peter.Thank you for your kind help.
----------
  I can't resist the temptation of presenting yet another variation on your integration problem, Zedong. It's an interesting problem. You will recall I derived a formula earlier whose result was called 'I' for the area integral of a polynomial term over a triangle (0,0), (x1,y1), (x2,y2) in terms of an I1, I2, and I3 which were for two triangles and a trapezoid and which together constituted the integral for just the desired triangle.

  Another way of deriving the same quantity, 'I', is to use iterated integration in which the inner integration takes place along rays emanating from the origin, with the result of this then subsequently integrated along the line segment connecting points (x1,y1) and (x2,y2). Without going into the details of this derivation, the result can be expressed as a single integral as follows. (We assume you are integrating the single term x^a*y^b for whole numbers a and b over the above triangle.)

 I = (x1*y2-x2*y1)/(a+b+2) * ...
     int('(x1+(x2-x1)*s)^a*(y1+(y2-y1)*s)^b','s',0,1);

The four arguments of 'int' are the integrand, the variable of integration s, and the lower and upper limits of s. The quantity s acts as a parameter which varies from 0 to 1 along the line segment. The quantities appearing outside the integral have been factored out because they remain constant as s varies throughout the integration. As you see. there is no inner integral here, the above integrand being already derived from such an inner "ray" integration process.

  If you use the Symbolic Toolbox to evaluate this, it will give a precise formula for the integral of f(x,y) = x^a*y^b over the triangle, just as my 'I' did in the previous article, and these should of course be mathematically identical expressions. The one here requires only one call on 'int'. However, the most important advantage of this approach in my opinion is that it is convenient for doing numerical integration in case you are dealing with large exponents a and b producing excessively long expressions from the 'int' function. In that case you could replace the above by a call to, say, 'quad', to evaluate it numerically. It is in a good form for numerical integration.

  Peter has shown you how to set up integration numerically based on doing the inner integration along the x-axis direction and the outer one in the y-direction, which he describes in terms of Green's Theorem. The one here differs in that the inner implicit integration is along rays from the origin and the outer integration along the far sides of the triangles. Both methods should give the same correct results, but I thought you would appreciate having a choice.

  Let me know if you wish to have the above integration formula, involving rays from the origin, explained.

Roger Stafford

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.

Contact us