Problem with quad2d integration

Question

meryem 2014-12-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/166422-problem-with-quad2d-integration

关闭： MATLAB Answer Bot 2021-8-20

Hi! I have a question on the quad2d function. When I run the code below ,it returns Q= 7.1754e-045. But it must be approximately 4.7399 e-043 . There is very little difference but it is very important for my project. When I decrease the value of 'AbsTol' ( for example 1e-85), it returns 4.867 e-043 but the code runs very slowly. I need another way to be fast and approximately 4.7399 e-043. Possible ?

fun = @(b,a) (b ./a).^19.*  exp(-(439 ./ a).^b) .* ((439 ./ a).^(b-1)) ...  .*  1./a .* 1./69.5391;
Q = quad2d(fun,6.00018,75.5392,248.063,1573.73 , 'AbsTol',1e-15,'Singular',false);

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Geoff Hayes 2014-12-12

Meryem - when I run the above code (R2014a, OS X 10.85) Q is 2.78617970352069e-21. Why is it necessary that your result be closer to 4.7399e-043? Please also elaborate on what you mean by the code runs very slowly - does this mean that it takes seconds, minutes, hours? Also, please indicate which version of MATLAB that you are using and the OS.

此问题已关闭。

Answer 1

Mike Hosea 2014-12-13

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/166422-problem-with-quad2d-integration#answer_162221

在 MATLAB Online 中打开

Is there a typo somewhere above? Here's a script I called "doit.m":

format short g
fun = @(b,a) (b ./a).^19.*  exp(-(439 ./ a).^b) .* ((439 ./ a).^(b-1)) .*  1./a .* 1./69.5391;
tic;
Qquad2dNS = quad2d(fun,6.00018,75.5392,248.063,1573.73,'RelTol',100*eps,'AbsTol',1e-50,'Singular',false,'MaxFunEvals',10000)
toc
tic;
Qquad2dS = quad2d(fun,6.00018,75.5392,248.063,1573.73,'RelTol',100*eps,'AbsTol',1e-50,'Singular',true,'MaxFunEvals',10000)
toc
tic;
QIterated = integral2(fun,6.00018,75.5392,248.063,1573.73,'RelTol',100*eps,'AbsTol',1e-50,'method','iterated')
toc
tic;
QTiled = integral2(fun,6.00018,75.5392,248.063,1573.73,'RelTol',100*eps,'AbsTol',1e-50,'method','tiled') 
toc
% Break the region up into n^2 pieces and integrate over each piece.
n = 30; 
x = linspace(6.00018,75.5392,n); 
y = linspace(248.063,1573.73,n); 
Q = 0; 
for i = 1:numel(x)-1
    for j = 1:numel(y)-1
        Q = Q + integral2(fun,x(i),x(i+1),y(j),y(j+1),'RelTol',100*eps,'AbsTol',1e-50,'method','iterated');
    end
end
QMIterated = Q
Q = 0; 
for i = 1:numel(x)-1
    for j = 1:numel(y)-1
        Q = Q + integral2(fun,x(i),x(i+1),y(j),y(j+1),'RelTol',100*eps,'AbsTol',1e-50,'method','tiled');
    end
end
QMTiled = Q
Answers = [QIterated,QTiled,QMIterated,QMTiled,Qquad2dNS,Qquad2dS]
Max = max(Answers);
Min = min(Answers);
ULPs = (Max - Min)/eps(Max)
>> doit
Qquad2dNS =
   2.0639e-18
Elapsed time is 0.763301 seconds.
Qquad2dS =
   2.0639e-18
Elapsed time is 1.608291 seconds.
QIterated =
   2.0639e-18
Elapsed time is 0.519476 seconds.
QTiled =
   2.0639e-18
Elapsed time is 1.408459 seconds.
QMIterated =
   2.0639e-18
QMTiled =
   2.0639e-18
Answers =
   2.0639e-18   2.0639e-18   2.0639e-18   2.0639e-18   2.0639e-18   2.0639e-18
ULPs =
    46

There are only 46 units of roundoff difference between the max and min for all these computed results, and the relative tolerance of 100*eps is only trying to get within 100 of the true solution. If that's not mathematically correct, then probably the numerical function isn't being evaluated accurately, which wouldn't surprise me given the large exponents.