Adaptive Simpson and Romberg integration methods

[robve]

My guess is that VA uses Guass Kronrod with pre-computed tables of weights and abscissas that can give at up to 100-decimal digit precision. The highlighted part I saw somewhere in his comments, so...

Gauss Kronrod for I1 gives -0.946030174332872 with 735 evaluations and \( 10^{-15} \) precision:

auto f2 = [](double x) { return cos(x)*log(x); }
double Q = gauss_kronrod<double, 15>::integrate(f2, 0, 1, 5, 1E-15, &error)

See NR2ed Ch.4.5 and Boost Math which I've used for this example.

PS: for clarity, 100 digits refers to the internal Gauss-Kronrod points, not the integral precision.

- Rob

Edit: fixed ugly "Gauss" typo.

[Valentin Albillo]

.
To:  robve and Albert Chan:

(03-26-2021 04:40 PM)robve Wrote:  I noticed that all of the integrals in his post have two distinguishing features: they have a U shape (or inverted) on both or either end that makes them suitable to apply a transform before numerical integration, i.e. change of variable or u-substitution. The trick to pull this off in an implementation is basically to have an initial check of a few points to find evidence of a U on either end and then transform

Seems to me we're talking different subject here.

On the one hand, I'm talking about good quadrature methods that are simply given the function to integrate, the limits of integration and the accuracy to achieve and upon being run the method returns the result meeting the specified accuracy as fast as possible. Period.

The user does nothing else other than supply the function, limits, and accuracy. He does not previously graph and analyze the function to see its shape, he does not study how to tame it if it needs taming, he does not do this or that transformation or another.

He just supplies the function, limits and accuracy and nothing else, The method is the one which does all the work and returns a reliable result, fast. That's in the spirit of what Mr. Kahan says in the HP Journal, I quote (my highlighting):

• "[...] such a procedure must be a computer program - call it P - that accepts as data two numerical values x and y and a program that calculates f(u) for any given value u, and from that data P must estimate I=∫(x,y),f(u) du. The integration procedure P is not allowed to read and understand the f-program but merely to execute it finitely often, as often as P likes, with any arguments u that P chooses."

My quadrature procedure P does exactly that, with the results seen in the 7 examples I posted.

On the other hand, you're doing exactly the opposite: you're using some quadrature methods which can't deal with difficult integrands and improper integrals in general, and you are doing the "read and understand"-ing of the f(u) on behalf on the quadrature method, graphing the curve, inspecting it for spikes or singularities, etc., then applying some transformation decided by you (not the quadrature method) to spoon-feed it to it, such that it can cope and deliver some acceptable results in acceptable times. Then you go on and say that the method is "really, really good" !?

Seriously: Are you joking ? Seems to me you're in a state of delusion, trying everything under the sun to reach the conclusion you want to reach (that the method is "really, really good"), instead of just throwing the function and the limits and desired accuracy to the method and letting it do what it can on its own, and if it cannot deliver reasonably without extensive help from the user, then you must simply recognize the fact and that's all, the world doesn't end.

Had I provided that kind of help to my procedure for I1, it wouldn't need 512 evaluations, just 16 at most, but I didn't do it because I'm trying to test the method, not my mathematical ability to help it cope. Matter of fact, the very first integral, I1, can be proved to be equivalent to Sin(x)/x, which is orders of magnitude easier to integrate and has only a removable singularity at x = 0, but I didn't apply any transformations and instead let my procedure cope with it on its own, which it did.


To wit, if this conversation is to be of any use, it's essential that we stick to a common ground, namely: you posit whatever methods you want to discuss but you let the methods do the work on their own, supplying them just the original f(x), not some transformation of it, the limits and the accuray, to afterwards judge them on their results and their speed.

To begin with, try my 7 examples as-is (with no transformations or tricks whatsoever) using your preferred methods and post the results here so that we can discuss them and compare their respective performances.

Agreed ?

Regards.
V.

[robve]

(03-25-2021 11:57 PM)Valentin Albillo Wrote:  Note: In everything that follows, no offence is ever meant, it's just my trademark style when I feel skeptic (let's say) about some statement or other, in particular if not supported by facts. IMHO, of course.

No offense. As I told my former students, critical thinking is important and do not stop questioning, but use the answers to analyze ideas and adjust them accordingly. We all build on existing work published by others and give credit for their efforts.

I am not sure what to think of statements like that and the references to "my method", which I never have used. I noticed that "your method" pzer on the HP-71B in every detail implements the standard Weierstrass / Durand-Kerner method for polynomial complex root finding aka "root polishing". I easily spotted this, because of the way you've set up the initial roots and the root polishing. Yet, nowhere is Durand-Kerner mentioned in your article. Further, it is the most bare and simple implementation.

There is actually more we can do beyond the basics, such as 1) removing trailing (in the code, leading in the poly) zero coefficients before polishing the roots, since these do not contribute to the result; 2) reject insignificant nonzero real or imaginary parts in the roots finally found. These require a little bit of code and incur almost no overhead. These additions are quite easy and often found in implementations.

So no offense as my post is derived from public existing work anyway, just like yours is.

- Rob

[robve]

(03-26-2021 05:52 PM)robve Wrote:  My guess is that VA uses Guass Kronrod with pre-computed tables of weights and abscissas that can give at up to 100-decimal digit precision. The highlighted part I saw somewhere in his comments, so...

Definitely just Gauss-Kronrod or/with double-exponential quadrature e.g. using tanh-sinh or similar:

The tanh_sinh integration routine will use this additional information internally when performing range transformation, so that for example, if one end of the range is zero (or infinite) then our transformations will get arbitrarily close to the endpoint without precision loss.

This gives the exact answer with only 74 evaluations: -0.946083070367183

auto f2 = [](double x) {return cos(x)*log(x); };
tanh_sinh<double> integrator;
double Q = integrator.integrate(f2, 0.0, 1.0);

The choice of tanh-sinh is usually, but not always, appropriate. The nice thing is that there are no bells and whistles, no parameters to figure out.

Again, this is not new and the method(s) are suitable for handling integrals with either singularities or large values near one or both endpoint(s), i.e. the U shape Albert commented on that we referred to earlier. Such shape issues can be resolved in several ways, sometimes by transformation or by picking a suitable integrator that does it (semi)-automatically.

PS. For info on double exponential integration and its implementation on the wp34S, see:
https://www.hpmuseum.org/forum/thread-8021.html

- Rob

Edit: added link to tanh-sinh and HP forum link.

[Valentin Albillo]

 
To: robve

(03-26-2021 06:41 PM)robve Wrote:  

(03-25-2021 11:57 PM)Valentin Albillo Wrote:  Note: In everything that follows, no offence is ever meant, it's just my trademark style when I feel skeptic (let's say) about some statement or other, in particular if not supported by facts. IMHO, of course.

No offense. [...]

I am not sure what to think of statements like that and the references to "my method", which I never have used. I noticed that "your method" pzer on the HP-71B in every detail implements the standard Weierstrass / Durand-Kerner method for polynomial complex root finding aka "root polishing". I easily spotted this, because of the way you've set up the initial roots and the root polishing. Yet, nowhere is Durand-Kerner mentioned in your article. Further, it is the most bare and simple implementation.

What does this have to do with the integration methods we were discussing in this thread  ??!!?

Is this intended as some kind of retaliatory ad-hominem attack by implying I don't give due credit or something ? Have you read my 40+ articles and noticed how they're being construed for the last 20+ years ?

Why would you do this, other than because you did take offence despite you denying it ?

For the benefit of the readers, robve sent me a PM almost a month ago about this and I kindly replied the following to him on 03-10-2021, 12:06 AM, that he seems to have forgotten. I quote my PM reply to him:

• "My articles are intended as just that, articles to be published in a physical fan-made magazine [...] or else on the WWW (MoHPC) for all kinds of fans, most of them not scholars, so my articles do not have the structure nor goals of a formal peer-reviewed paper.
  
  This means I usually include a brief Introduction, then the program Code, usage Instructions, several fully-worked Examples an some Notes.
  
  This already takes vast amounts of my free time (about 1 month per article, from the idea to the finished PDF ready for publication, including testing) so I simply can't allocate more time to also include Algorithms, Formulae and References, as I would do for a formal paper intended for a peer-reviewed publication and a vastly different target audience."

Seems robve forgot or simply didn't want to ruin his uncalled-for ad-hominem attack by including my reply here.

Quote:There is actually more we can do beyond the basics, such as 1) removing trailing (in the code, leading in the poly) zero coefficients before polishing the roots, since these do not contribute to the result; 2) reject insignificant nonzero real or imaginary parts in the roots finally found. These require a little bit of code and incur almost no overhead. These additions are quite easy and often found in implementations.

Again, this is but a very poorly disguised attack on the quality of my implementations by mentioning that they're basic and can be easily improved. He seems to have conveniently forgotten (again!) my reply to this in the very same PM or else he doesn't see fit to include it here, so I'll do it myself. I replied to him:

• "One other thing: as creating an article takes so much time, once created I don't allocate any more resources to try and improve it or research it in any way, I simply jump to the next article because that's the only way I can get things moving forward, else I would still be researching and optimizing my article Minimax Polynomial Fit (recommended).
  
  Right now I've completely put aside the PZER 71B article and I'm working on two other unrelated ones simultaneously, which right now are 50% and 25% complete so they'll be uploaded and announced next April/May.
  
  In short: I don't optimize or allocate more time to past articles. I dedicate a fair amount of time to them while creating them but afterwards (except to correct blatant errors) I have no time for them at all, so it's really fortunate that you're able to consider optimizations and I thank you for your interest in experimenting with them. Actually that's the whole idea: I publish something and interested people have a fair chance at non-trivial optimizations, having a good time while at it and possibly improving their knowledge and skills."

Also, I see that you haven't read or at any rate don't comment on my on-topic (integration !!) post #22, but that's irrelevant now.

Listen, I'm no fool and can clearly see by this post of yours a very many things that I'm itching to say but can't and won't discuss in public lest they'd be considered ad-hominens, but also I'm not a person to suffer them gladly as well, nor the kind of individuals who issue them, so please take note of the following:

• I'm out of this thread. Whatever you may or may not reply, I won't read.
   
• Please do not PM me in the future. I won't reply to or reference any posts or PMs from you, except to notify the moderators in case you disrespect me.
   
• I was completing a new Short & Sweet Math Challenge #26 to be published within a week and two new Articles to be published on April 15 and May 1, respectively, but having suffered this more than once in the past, I am so fed up with people like you that I won't publish them and further I will quit for a while from the forum till the disgust recedes at seeing all my disinterested hard work and well-meant efforts being questioned with blatant bad intent by the ones like you.
  

V.

[Gene]

I am calling a time out on this thread.

Gene