How to calculate a double integral inside the domain of intersection of two functions?

I want to calculate a double integral of an arbitrary function (f) inside the region of intersection of two other functions. Please suggest a fast and convenient approach.

clear
JJ = 5;
II = 5;
W = rand(II, JJ);
syms x y
w = sym('0');
f = sym('0');
for i=1:II
    for j=1:JJ
        w =w+W(i, j)*legendreP(i-1, x)*legendreP(j-1, y);
        f =f+(legendreP(i-1, x)*legendreP(j-1, y))^2;
    end
end
H = 0.5*(1+tanh(w));
fsurf(w,[-1,1,-1,1],'red')
hold on
% figure
fsurf(H,[-1,1,-1,1],'blue')
%F = double integral of f inside the domain of intersection of two functions as
%the region showed in pic

[Matt J]

Is the area shaded in black the region to be integrated? It is not clear what defines that region. The red and black surfaces do not seem to be touching there, so in what sense do they "intersect"?

[Mehdi]

Yes the black shaded region is the domain that the integral must be computed. As I said this region is the upper region created from the intersection of two functions H and w in -1<x,y<1.

[Torsten]

You mean the region somehow enclosed by the curve H-w = 0 ?

Then integrate f.*(H<w) over [-1,1] x [-1,1].

rng("default")
JJ = 5;
II = 5;
W = rand(II, JJ);
syms x y
w = sym('0');
for i=1:II
    for j=1:JJ
        w =w+W(i, j)*legendreP(i-1, x)*legendreP(j-1, y);
    end
end
H = 0.5*(1+tanh(w));
Hminusw = fimplicit(matlabFunction(H-w),[-2 2 -2 2]);

[Mehdi]

Yes, the double integral of f inside the region enclosed by the curve H-w = 0 and (H<w) over [-1,1] x [-1,1], as I showed by red color in following figure.

[Torsten]

Then do as I suggested:

Define D = f*(H<w) and integrate D over [-1,1]^2.

Note that H<w gives 1 if H<w and 0 else. Thus you integrate f over the area in red.

[Mehdi]

What do mean by H<w gives 1 if H<w and 0 else? Can you show me on the codes?

[Torsten]

x = [1 2 3 4];
y = [3 4 1 2];
x<y

ans = 1×4 logical array
   1   1   0   0

% Compute area of circle with radius 1
f = @(x,y) (x.^2+y.^2<=1)*1;
value = integral2(f,-2,2,-2,2)

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

value = 3.1415

[Mehdi]

Could not get answer!

clear
JJ = 5;
II = 5;
W = rand(II, JJ);
syms x y
w = sym('0');
f = sym('0');
for i=1:II
    for j=1:JJ
        w =w+W(i, j)*legendreP(i-1, x)*legendreP(j-1, y);
        f =f+(legendreP(i-1, x)*legendreP(j-1, y))^2;
    end
end
H = 0.5*(1+tanh(w));
D = simplify(f)*(H<w) ;
Dint= int(int(D ,x,-1,1),y,-1,1)

Warning: Unable to check whether the integrand exists everywhere on the integration interval.

Dint = 

NaN

fint=int(int(f ,x,-1,1),y,-1,1)

fint = 

[Torsten]

I don't know how trustworthy the results are, but I don't see a different way to get what you want.

But the results don't seem that bad since added, they give the integral of f over [-1,1]^2.

rng("default")
JJ = 5;
II = 5;
W = rand(II, JJ);
syms x y
w = sym('0');
f = sym('0');
for i=1:II
    for j=1:JJ
        w =w+W(i, j)*legendreP(i-1, x)*legendreP(j-1, y);
        f =f+(legendreP(i-1, x)*legendreP(j-1, y))^2;
    end
end
H = 0.5*(1+tanh(w));
f = matlabFunction(f)

f = function_handle with value:
    @(x,y)y.^2.*(x.*(3.0./2.0)-x.^3.*(5.0./2.0)).^2+x.^2.*(y.*(3.0./2.0)-y.^3.*(5.0./2.0)).^2+x.^2.*y.^2+(x.*(3.0./2.0)-x.^3.*(5.0./2.0)).^2.*(y.^2.*(3.0./2.0)-1.0./2.0).^2+(y.*(3.0./2.0)-y.^3.*(5.0./2.0)).^2.*(x.^2.*(3.0./2.0)-1.0./2.0).^2+(x.^2.*(-1.5e+1./4.0)+x.^4.*(3.5e+1./8.0)+3.0./8.0).^2+x.^2.*(y.^2.*(3.0./2.0)-1.0./2.0).^2+y.^2.*(x.^2.*(3.0./2.0)-1.0./2.0).^2+(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0).^2+(x.^2.*(3.0./2.0)-1.0./2.0).^2.*(y.^2.*(3.0./2.0)-1.0./2.0).^2+(x.*(3.0./2.0)-x.^3.*(5.0./2.0)).^2.*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0).^2+(y.*(3.0./2.0)-y.^3.*(5.0./2.0)).^2.*(x.^2.*(-1.5e+1./4.0)+x.^4.*(3.5e+1./8.0)+3.0./8.0).^2+y.^2.*(x.^2.*(-1.5e+1./4.0)+x.^4.*(3.5e+1./8.0)+3.0./8.0).^2+x.^2.*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0).^2+(y.^2.*(3.0./2.0)-1.0./2.0).^2.*(x.^2.*(-1.5e+1./4.0)+x.^4.*(3.5e+1./8.0)+3.0./8.0).^2+(x.^2.*(3.0./2.0)-1.0./2.0).^2.*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0).^2+(x.^2.*(-1.5e+1./4.0)+x.^4.*(3.5e+1./8.0)+3.0./8.0).^2.*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0).^2+(x.*(3.0./2.0)-x.^3.*(5.0./2.0)).^2+(y.*(3.0./2.0)-y.^3.*(5.0./2.0)).^2+x.^2+y.^2+(x.^2.*(3.0./2.0)-1.0./2.0).^2+(y.^2.*(3.0./2.0)-1.0./2.0).^2+(x.*(3.0./2.0)-x.^3.*(5.0./2.0)).^2.*(y.*(3.0./2.0)-y.^3.*(5.0./2.0)).^2+1.0

g = matlabFunction(H-w)

g = function_handle with value:
    @(x,y)x.*4.642718471329099e-1+y.*1.152891029414135e-1+tanh(x.*(-4.642718471329099e-1)-y.*1.152891029414135e-1+y.*(x.^2.*8.203222788074758e-1-2.734407596024919e-1)-y.*(x.*1.436260253151446-x.^3.*2.393767088585744)-x.*(y.*(3.0./2.0)-y.^3.*(5.0./2.0)).*4.21761282626275e-1+x.*y.*2.784982188670484e-1-(y.*(3.0./2.0)-y.^3.*(5.0./2.0)).*(x.^2.*(-3.598096598973386)+x.^4.*4.197779365468951+3.598096598973386e-1)+(y.^2.*(3.0./2.0)-1.0./2.0).*(x.^2.*(-3.001051758333)+x.^4.*3.5012270513885+3.001051758333e-1)+x.*(y.^2.*(3.0./2.0)-1.0./2.0).*9.705927817606157e-1+y.*(x.^2.*(-3.618332006997287)+x.^4.*4.221387341496835+3.618332006997287e-1)-(x.*7.280634730842618e-1-x.^3.*1.213439121807103).*(y.^2.*(3.0./2.0)-1.0./2.0)-(x.*1.400989871636326-x.^3.*2.334983119393876).*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0)+x.*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0).*3.571167857418955e-2+(x.^2.*1.273693958803166-4.245646529343886e-1).*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0)+(x.^2.*(-2.545256830716651)+x.^4.*2.969466302502759+2.545256830716651e-1).*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0)-x.^2.*2.180866948905027+x.^3.*2.283439640347548+x.^4.*2.766571702236167-y.^2.*2.222607999320878+y.^3.*3.547158465680383e-1+y.^4.*2.868865558810067+(y.^2.*(3.0./2.0)-1.0./2.0).*(x.^2.*1.435750422364418-4.785834741214728e-1)-(x.^2.*1.373603287783601-4.578677625945335e-1).*(y.*(3.0./2.0)-y.^3.*(5.0./2.0))+(x.*1.188310994339332-x.^3.*1.980518323898886).*(y.*(3.0./2.0)-y.^3.*(5.0./2.0))+1.1554612169259)./2.0-y.*(x.^2.*8.203222788074758e-1-2.734407596024919e-1)+y.*(x.*1.436260253151446-x.^3.*2.393767088585744)+x.*(y.*(3.0./2.0)-y.^3.*(5.0./2.0)).*4.21761282626275e-1-x.*y.*2.784982188670484e-1+(y.*(3.0./2.0)-y.^3.*(5.0./2.0)).*(x.^2.*(-3.598096598973386)+x.^4.*4.197779365468951+3.598096598973386e-1)-(y.^2.*(3.0./2.0)-1.0./2.0).*(x.^2.*(-3.001051758333)+x.^4.*3.5012270513885+3.001051758333e-1)-x.*(y.^2.*(3.0./2.0)-1.0./2.0).*9.705927817606157e-1-y.*(x.^2.*(-3.618332006997287)+x.^4.*4.221387341496835+3.618332006997287e-1)+(x.*7.280634730842618e-1-x.^3.*1.213439121807103).*(y.^2.*(3.0./2.0)-1.0./2.0)+(x.*1.400989871636326-x.^3.*2.334983119393876).*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0)-x.*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0).*3.571167857418955e-2-(x.^2.*1.273693958803166-4.245646529343886e-1).*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0)-(x.^2.*(-2.545256830716651)+x.^4.*2.969466302502759+2.545256830716651e-1).*(y.^2.*(-1.5e+1./4.0)+y.^4.*(3.5e+1./8.0)+3.0./8.0)+x.^2.*2.180866948905027-x.^3.*2.283439640347548-x.^4.*2.766571702236167+y.^2.*2.222607999320878-y.^3.*3.547158465680383e-1-y.^4.*2.868865558810067-(y.^2.*(3.0./2.0)-1.0./2.0).*(x.^2.*1.435750422364418-4.785834741214728e-1)+(x.^2.*1.373603287783601-4.578677625945335e-1).*(y.*(3.0./2.0)-y.^3.*(5.0./2.0))-(x.*1.188310994339332-x.^3.*1.980518323898886).*(y.*(3.0./2.0)-y.^3.*(5.0./2.0))-6.554612169259004e-1

D1 = @(x,y)f(x,y).*(g(x,y)<0);
D2 = @(x,y)f(x,y).*(g(x,y)>0);
D1int = integral2(D1,-1,1,-1,1)

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

D1int = 5.7859

D2int = integral2(D2,-1,1,-1,1)

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

D2int = 6.9919

D1int + D2int

ans = 12.7778

Fint = integral2(f,-1,1,-1,1)

Fint = 12.7778

[Mehdi]

nice!

In the case which w becomes zero, D faced error. why?

clear
JJ = 5;
II = 5;
W = 0*rand(II, JJ);
syms x y
w = sym('0');
f = sym('0');
for i=1:II
    for j=1:JJ
        w =w+W(i, j)*legendreP(i-1, x)*legendreP(j-1, y);
        f =f+(legendreP(i-1, x)*legendreP(j-1, y))^2;
    end
end
H = 0.5*(1+tanh(w));
f = matlabFunction(f);
g = matlabFunction(H-w);
D = @(x,y)f(x,y).*(g(x,y)<0);
Dint = integral2(D,-1,1,-1,1)

Error using symengine>@()1.0./2.0
Too many input arguments.

Error in solution (line 18)
D = @(x,y)f(x,y).*(g(x,y)<0);

Error in integral2Calc>integral2t/tensor (line 228)
Z = FUN(X,Y); NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Fint = integral2(f,-1,1,-1,1)

[Torsten]

Because g is constant and MatlabFunction turns g into a numerical function with no input arguments.

You might use

JJ = 5;
II = 5;
W = 0*rand(II, JJ);
syms x y
w = sym('0');
f = sym('0');
for i=1:II
    for j=1:JJ
        w =w+W(i, j)*legendreP(i-1, x)*legendreP(j-1, y);
        f =f+(legendreP(i-1, x)*legendreP(j-1, y))^2;
    end
end
H = 0.5*(1+tanh(w));
f = matlabFunction(f,'Vars',[x y]);
g = matlabFunction(H-w,'Vars',[x y]);
D = @(x,y)f(x,y).*(g(x,y)<0);
Dint = integral2(D,-1,1,-1,1)
Fint = integral2(f,-1,1,-1,1)

instead.

[Mehdi]

Any suggestion to speed up the calculation?

clc
clear
JJ = 9;
II = 9;
W = 01*rand(II, JJ);
syms x y
w = sym('0');
f = sym('0');
for i=1:II
    for j=1:JJ
        w =w+W(i, j)*legendreP(i-1, x)*legendreP(j-1, y);
        f =f+(legendreP(i-1, x)*legendreP(j-1, y))^2;
    end
end
H = 0.5*(1+tanh(11*w));
f = matlabFunction(f,'Vars',[x y]);
g = matlabFunction(H-0.5,'Vars',[x y]);
D = @(x,y)f(x,y).*(g(x,y)>0);
Dint = integral2(D,-1,1,-1,1)

[Torsten]

Any suggestion to speed up the calculation?

No. I remember we had this integration problem before. And the region where g(x,y)>0 is quite complicated.

Answer 1

[Bjorn Gustavsson]

Perhaps you can use Green's theorem (you'd be very lucky if you could - but if you were to be that lucky in this case it would be a shame to miss it). That would take you from a sum of integrals over rather complicated regions to perhaps simpler integrals around the boundaries of the region. That would be nice. Given the shape of your function it doesn't seem entirely improbable.

For the case where you actually have to perform the calculations you would use the steps suggested in @Torsten's comment.

HTH

[Mehdi]

I think it is impossible to find M and L to use Green's theorem.

[Torsten]

M and L might work - the main problem is C.

[Bjorn Gustavsson]

@Torsten: I thought(hoped) it would be possible to piece together C from the output of fimplicit (and the [0,1], [0,1] edges). Then it ought to be a "reasonably" straightforward number of integrals.

[Mehdi]

C can be found by FEX submission availiable through https://www.mathworks.com/matlabcentral/fileexchange/74010-getcontourlinecoordinates?s_tid=srchtitle. The problem for me is to find M and L (a,b). Any suggestions?

[Torsten]

a = b = 0.

You can use any function M with dM/dx = f.

But you need a symbolic integration here.

I doubt the whole process will be faster than using integral2 directly, apart from the problems of evaluating the curve integral.