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

Question

Mehdi 2023-1-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1895195-how-to-calculate-a-double-integral-inside-the-domain-of-intersection-of-two-functions

编辑： Torsten 2023-1-25

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

13 个评论
显示 11更早的评论隐藏 11更早的评论

Torsten 2023-1-17

编辑：Torsten 2023-1-17

在 MATLAB Online 中打开

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 2023-1-18

在 MATLAB Online 中打开

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 2023-1-18

编辑：Torsten 2023-1-18

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 2023-1-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1895195-how-to-calculate-a-double-integral-inside-the-domain-of-intersection-of-two-functions#answer_1150380

编辑：Bjorn Gustavsson 2023-1-17

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