This post is the answer to this question.
The procedure named IntOverDomain finds a double integral over an arbitrary domain bounded by a non-selfintersecting piecewise smooth curve. The code of the procedure uses the well-known Green's theorem.
Each section in the border should be specified by a list in the following formats :
1. If a section is given parametrically, then [[f(t), g(t)], t=t1..t2]
2. If several consecutive sections of the border or the entire border is a broken line, then it is sufficient to set vertices of this broken line [ [x1,y1], [x2,y2], .., [xn,yn] ] (for the entire border should be [xn,yn]=[x1,y1] ).
Required parameters of the procedure: f is an expression in variables x and y , L is the list of all the sections. The sublists of the list L must follow in the positive direction (counterclockwise).
The code of the procedure:
restart;
IntOverDomain := proc(f, L)
local n, i, j, m, yk, yb, xk, xb, Q, p, P, var;
n:=nops(L);
Q:=int(f,x);
for i from 1 to n do
if type(L[i], listlist(algebraic)) then
m:=nops(L[i]);
for j from 1 to m-1 do
yk:=L[i,j+1,2]-L[i,j,2]; yb:=L[i,j,2];
xk:=L[i,j+1,1]-L[i,j,1]; xb:=L[i,j,1];
p[j]:=int(eval(Q*yk,[y=yk*t+yb,x=xk*t+xb]),t=0..1);
od;
P[i]:=add(p[j],j=1..m-1) else
var := lhs(L[i, 2]);
P[i]:=int(eval(Q*diff(L[i,1,2],var),[x=L[i,1,1],y=L[i,1,2]]),L[i,2]) fi;
od;
add(P[i], i = 1 .. n);
end proc:
Examples of use.
1. In the first example, we integrate over a quadrilateral:
with(plottools): with(plots):
f:=x^2+y^2:
display(polygon([[0,0],[3,0],[0,3],[1,1]], color="LightBlue")); # Visualization of the domain of integration
IntOverDomain(x^2+y^2, [[[0,0],[3,0],[0,3],[1,1],[0,0]]]); # The value of integral
2. In the second example, some sections of the boundary of the domain are curved lines:
display(inequal({{y<=sqrt(x),y>=sin(Pi*x/3)/2,y<=3-x}, {y>=-2*x+3,y>=sqrt(x),y<=3-x}}, x=0..3,y=0..3, color="LightGreen", nolines), plot([[t,sqrt(t),t=0..1],[t,-2*t+3,t=0..1],[t,3-t,t=0..3],[t,sin(Pi*t/3)/2,t=0..3]], color=black, thickness=2));
f:=x^2+y^2: L:=[[[t,sin(Pi*t/3)/2],t=0..3],[[3,0],[0,3],[1,1]], [[t,sqrt(t)],t=1..0]]:
IntOverDomain(f, L);
3. If f=1 then the procedure returns the area of the domain:
IntOverDomain(1, L); # The area of the above domain
evalf(%);
IntOverDomain.mw
Edit.