How can i know the size of the rectangles in a numerical integration?

4 views (last 30 days)
Pau
Pau on 5 Nov 2013
Edited: Mike Hosea on 6 Nov 2013
I am doing some simulations in optics and i need to know how the function that evaluates my integral discretize the area of integration. I am looking for this in the code of quad2d and integral2, but i can not find it.
Thanks

Answers (2)

Mike Hosea
Mike Hosea on 5 Nov 2013
Edited: Mike Hosea on 5 Nov 2013
All the INTEGRAL and QUAD* functions are adaptive quadrature routines. INTEGRAL2 with the 'tiled' method transforms the integration region to a rectangle, including a transform that weakens edge singularities. In the transformed domain the rectangle is subdivided into 4 equal rectangles, each of which may be subdivided into 4 equal rectangles if the error estimate for that part indicates the need, and so forth. It uses a tensor product Gauss-Kronrod rule on each rectangle to approximate the integral over that rectangle and to estimate the error. QUAD2D is the same as INTEGRAL2, but QUAD2D also provides an option not to include the singularity-weakening component of the transform. You can read about the algorithm in the reference (an article by L.F. Shampine) found on the doc page of QUAD2D.
The 'iterated' method of INTEGRAL2 takes a different approach. As the name suggests, it evaluates the double integral as an iterated integral. Both inner and outer integrations are performed with the algorithm of INTEGRAL. INTEGRAL uses the same basic algorithm as QUADGK, which likewise has a reference on the doc page.
  3 Comments
Mike Hosea
Mike Hosea on 6 Nov 2013
To answer the question posed below, the "FailurePlot" plots only the pieces that still need refining because the error estimate is not small enough over those areas. The pieces that have passed the error test do not appear. What you see below would be the early stages of adaptive quadrature, so initially, yes, that is what it is doing. The final division, however, will probably not be so uniform, but you can visualize it as sub-dividing some of those pieces into 4 new pieces and not others, and then sub-dividing some of those new ones into 4 new pieces and not other, etc. QUAD2D does this with the goal of driving down approximate upper bound on the global error until it satisfies errest <= max(abstol,reltol*abs(Q)), where Q is the current approximation of the integral.
PS: The transform itself introduces a mild singularity because circle is not mapped smoothly to a rectangle, so if you were really trying to solve a problem there, you would probably not want to set 'Singular' to false.

Sign in to comment.


Pau
Pau on 6 Nov 2013
Then if I get
using
ymin = @(x) -sqrt(a^2 - x.^2);
ymax = @(x) +sqrt(a^2 - x.^2);
P(i,1)=quad2d(fun,a,a,ymin,ymax,'AbsTol',10^-17,'MaxFunEvals',22,'FailurePlot',true,
'Singular',false);
This means that quad2d is subdividing my circle in this pieces?

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!