Perheps not a bug, but a… not properly working method of adaptive integral for my problem.
I atttach the code - the problem is the way of integrating in functions re_int and im_int. When I return result of manual integrating (intg variable calculation) it’s fine. When I use TF2::Integral, it’s wrong. Below is the description of the code, for it is difficult to make it any shorter.
I am calculating the image of a light source illuminating circular aperture and than passing through a thing lens. I use fresnel integral. In interiorre() and interiorim() I calculate the sine and cosine of the exp(-ikr) term. In re_int() and im_int() I calculate integrals of those functions, in psf_calc() I get the intensity from those integrals.
The result should be a maximum in the centre of the image - it is so for the manual integration. For the TF2::Integral, I get minimum in the centre of the image. psf2d.cpp (4.36 KB)
I have modified your code such that
-it compiles
-it runs in a reasonable amount of time (your own loops take for ever)
-it produces what I believe are reasonable results, but we cannot provide further assistance on algorithmic problems
Indeed, during a code cleanup I commented one thing to much and the code didn’t compile. I am sorry to waste your time with this. However, that is not the problem. The length of running is a normal thing.
I compare here results of ROOT integral: TF2:Integral() and a manual integration of TF2 in a loop (starting on line 46 and 70). Results are different - manual integration produces maximum in the centre of the canvas, while ROOT integral produces minimum. As far as I know, the result of the manual integration is the good one (I compare it to results of other programs).
I attach a new code. Please run the to options and compare:
.x psf2d.cpp+(false)
.x psf2d.cpp+(true)
the first one is the ROOT integration of TF2, the second (which here takes about a minute to calculate) is a manual loop calculation of the integral.
This may be due to numerical bug or simply uncertainty of ROOT adaptive integral algorithm… psf2d.cpp (3.69 KB)
Hi. I’m not versed at Geometric Optics but after seeing the imint and reint functions surface’s being drawn I think that perhaps both functions are too much of oscillating ones for the standard available numeric approximations.
A question that I have is since you consider fresnel integrals: the closest formula to the interiorre and interiorim ones that I have found takes the squared root of the sine (and cosine) argument(s) as denominator, not just the argument itself. Is this correct?
[quote=“brabacal”]Hi. I’m not versed at Geometric Optics but after seeing the imint and reint functions surface’s being drawn I think that perhaps both functions are too much of oscillating ones for the standard available numeric approximations.
[/quote]
It is possible. If this is true, I guess it would be nice to have option to make standard, slow integral for TF2.
If I understand you correctly - yes. However my interiorre and interiorim are still a rough attempt, failed on the integral problem…
Yes, we have been looking (in particular Bartolomeu) into this.
The numerical integration cannot handle correctly a function oscillating so much. You acns ee this by playing with the parameters of root.cern.ch/root/htmldoc/TF1.ht … alMultiple
To get a decent error you will need a very large number of points, which will require a large CPU time.
I would work on the Math of your Fresnel integrals. Normally you should be able to perform the integration in the complex planes and probably also express the result in terms of some known special functions, which exist in ROOT or other numerical libraries.
We don’t need to add a new method for this, since it is not very precise for a normal function and it will depend on the number of fixed sampling points.
You can do this easily by going through an histogram. Here is an example:
TF2 * f2 = new TF2("f2","x+y",0,1,0,1);
f2->Draw(); // needed to create the underlying histogram
double integral = f2->GetHistogram()->Integral("width");
and with TF2::SetNpx and TF2::SetNpy you can control the number of sampling points.
I can see your point. However using the TF2::Draw() is time consuming - manual integration using TF2::Eval() and suming up results is much faster here, although still slow. And I need to use integration results for fitting, so any speed improvement is important.
Anyway thanks for the whole discussion. The only thing I suggest is adding a remoark in the TF2::Integral() documentation, that this integration method may not work for fast oscilating functions. Could save time for numerical ignorants like me
It is true, doing directly could be a bit faster, although the I think the majority of time will be in evaluating the function value for creating the grid and that will be the same.
Thank you for the feedback, I will probably add a comment in the Integral method. I know there exist some numerical method for oscillating functions in 1D, maybe it will be worth integrating in ROOT.
I will investigate if exist something in multi-dimensions
