Path: news.mathworks.com!not-for-mail From: "zedong " <zdongwu@gmail.com> Newsgroups: comp.soft-sys.matlab Subject: Re: How to integration a polynomial on a polygon? Date: Mon, 22 Dec 2008 01:01:53 +0000 (UTC) Organization: The MathWorks, Inc. Lines: 55 Message-ID: <gimou1$q54$1@fred.mathworks.com> References: <gikgkj$4qp$1@fred.mathworks.com> <gil3rj$ss6$1@fred.mathworks.com> <gimie6$gr7$1@fred.mathworks.com> Reply-To: "zedong " <zdongwu@gmail.com> NNTP-Posting-Host: webapp-05-blr.mathworks.com Content-Type: text/plain; charset="ISO-8859-1" Content-Transfer-Encoding: 8bit X-Trace: fred.mathworks.com 1229907713 26788 172.30.248.35 (22 Dec 2008 01:01:53 GMT) X-Complaints-To: news@mathworks.com NNTP-Posting-Date: Mon, 22 Dec 2008 01:01:53 +0000 (UTC) X-Newsreader: MATLAB Central Newsreader 1639198 Xref: news.mathworks.com comp.soft-sys.matlab:508263 "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.