Dear Alexei, Thank you for your interest and careful examination of the code. I made an effort to fit the functions in the GiNaC framework and got series expansions etc. more or less for free. I'm looking forward to seeing the result when it meets the GiNaC standards. To answer your questions a fundamental property of the algorithm is needed: Depending on the argument a so called "steepest descent" contour is determined along which the integrand is monotonically decreasing. This ensures that a) roundoff does not occur, b) the range over which the integrand has to be evaluated to obtain a certain accuracy is minimized and c) the dominant factor (possibly very larger or small) can be placed outside the integral. Furthermore all derivatives at +/- infinity are zero. For functions with this property the trapezoidal rule has exponential convergence, that is, halving the stepsize more than doubles the number of significant digits of the result. For an integrator with non-constant step size this property would be lost.
AiryContour* pcntr=createAiryContour(zn);
Pointers are evil. What is the reason for explicit
allocation/deallocation?
Again depending on the argument (zn) but more involved, 3 different types of contours are needed of which one is the truly steepest descent contour and the other 2 are approximations thereof. AiryContour is the base class of them. I hope this answers your questions Best regards, Jan
Hello, On Sun, Dec 23, 2007 at 02:53:55PM +0100, Jan Bos wrote:
To answer your questions a fundamental property of the algorithm is needed:
Depending on the argument a so called "steepest descent" contour is determined along which the integrand is monotonically decreasing.
This ensures that a) roundoff does not occur,
Roundoff errors occur (almost) always, unless you've got a magic hardware which operates on real numbers (as opposed to floating point ones). However, in a carefully designed code these errors are smaller (or of the same order) than the errors of the algorithm itself (i.e. due to truncation of a series, replacing an integral with a sum, etc). I don't quite understand how your code solves this problem.
b) the range over which the integrand has to be evaluated to obtain a certain accuracy is minimized and
Anyway, the higher accuracy, the bigger this interval is. With IEEE floating point numbers (or any other numbers with the fixed size of the mantissa) arbitrary high accuracy does not make sense, so this is (probably) not a problem. However, for arbitrary precision calculation the accuracy can be arbitrary high, so the interval could be arbitrary big.
c) the dominant factor (possibly very larger or small) can be placed outside the integral.
Furthermore all derivatives at +/- infinity are zero. For functions with this property the trapezoidal rule has exponential convergence, that is, halving the stepsize more than doubles the number of significant digits of the result.
That would be correct in real number arithmetics. With floating point numbers halving the step size can result in bigger roundoff error (becouse of doubled number of operations). Thus, one need to find an optimal step size, so both roundoff error and the error of algorithm itself are small enough. If the target accuracy is very high[1], such a step size might not exist at all. Hence my question -- is numerical integration suitable for *arbitary* precision calculation? [1] and/or the integrand is weird, i.e. sin(1/x) when x -> 0
For an integrator with non-constant step size this property would be lost.
I don't quite understand. Actually, if the function is monotoneous along the integration contour, I'd expect adaptive integration algorithm to converge faster (and give smaller roundoff errors). So, could you elaborate, please? Best regards, Alexei -- All science is either physics or stamp collecting.
participants (2)
-
Alexei Sheplyakov
-
Jan Bos