Post Reply 
Numerical integration methods
07-25-2022, 10:36 PM
Post: #21
RE: Numerical integration methods
That is what is going on! When, on the Classpad 400, I write integral from 1 to 100 of x^(2(e^((1/x)^3)-1)), I get the 99.5996+ answer!
Find all posts by this user
Quote this message in a reply
07-26-2022, 05:18 PM
Post: #22
RE: Numerical integration methods
I never considered that as a possible interpretation of the original expression. Well done for figuring out what was going on! I'm glad that's sorted.

Nigel (UK)
Find all posts by this user
Quote this message in a reply
07-31-2022, 04:59 PM
Post: #23
RE: Numerical integration methods
(07-24-2022 08:55 PM)Wes Loewer Wrote:  (Note: It's possible that I am mistranslating "ordre". What the paper called "ordre 30" seems to be "degree 29" (ie, 30 terms), but I could be wrong about this. Francophones, please correct me if I am mistaken.)
"order" may mean max degree of polynomial where the method is exact, or the exponent in the step of the error, up to the author choice...

Quote:Since the 15-node Gauss-Kronrod is exact up to degree 22 and the Prime's method is exact up to degree 29, it would seem that the Prime's method is a bit better with only a small added overhead. A slightly more accurate calculation per iteration could occasionally reduce the need for further recursion. It would be interesting to do some comparisons of these two methods with some well-behaved functions as well as more "temperamental" ones.

On a final note, I am curious why Hairer chose only 6 nodes for the 3rd calculation instead of 7. Any insight on that would be welcomed. [Edit: Perhaps it is important to have the 3rd calculation's 6 nodes be a subset of the 2nd calculation's 14 nodes. This would not be the case if 7 nodes were used for the 3rd calculation. It is not clear to me why this subset would be necessary or even desirable.]
Yes, it is important not to re-evaluate the function for efficiency reasons. The idea is to estimate the error of the order 30 method in h^30 by
err1*(err1/err2)^2
where err1=abs(i30-i14) is in h^14 and err2=abs(i30-i6) in h^6
Find all posts by this user
Quote this message in a reply
08-01-2022, 03:43 AM
Post: #24
RE: Numerical integration methods
(07-31-2022 04:59 PM)parisse Wrote:  "order" may mean max degree of polynomial where the method is exact, or the exponent in the step of the error, up to the author choice...

I see. The order of the error, not the order of the polynomial, in this case. Thank you.

(07-31-2022 04:59 PM)parisse Wrote:  Yes, it is important not to re-evaluate the function for efficiency reasons. The idea is to estimate the error of the order 30 method in h^30 by
err1*(err1/err2)^2
where err1=abs(i30-i14) is in h^14 and err2=abs(i30-i6) in h^6

This method works beautifully, but I am fuzzy about some of the details.

For instance, why *(err1/err2)^2 ? Why not *(err1/err2) or *(err1/err2)^3 ? Was ^2 derived mathematically, or experimentally? I'm guessing it was experimentally determined to give a reasonable error approximation.

And why does the 3rd calculation use 6 nodes? Why not 5 or 7? I realize you want the 3rd calculation to be different enough (but not too different) from the 2nd calculation to determine the error. Was 6 determined experimentally to produce an optimal error estimate? Or is there something mathematical that makes 6 the best choice? Using 6 leaves the entire middle third of the interval unevaluated. Seems like including the center node would have been advantageous.
Find all posts by this user
Quote this message in a reply
08-01-2022, 03:01 PM
Post: #25
RE: Numerical integration methods
(08-01-2022 03:43 AM)Wes Loewer Wrote:  For instance, why *(err1/err2)^2 ? Why not *(err1/err2) or *(err1/err2)^3 ? Was ^2 derived mathematically, or experimentally? I'm guessing it was experimentally determined to give a reasonable error approximation.

And why does the 3rd calculation use 6 nodes? Why not 5 or 7?

This is my guess ...

err1 = k1 * h^15
err2 = k2 * h^7

integral error estimate = err1*(err1/err2)^2 = k1*(k1/k2)^2 * h^(15+8*2) = k * h^31

If (k, k1, k2) are similar in size, constant term matched too.

Exponents cannot be picked in random; it had to produce O(h^31) on the right.
If we use 5 or 7 nodes, instead of square, it would need some non-integer exponents.
Find all posts by this user
Quote this message in a reply
08-01-2022, 03:46 PM
Post: #26
RE: Numerical integration methods
(08-01-2022 03:01 PM)Albert Chan Wrote:  This is my guess ...

err1 = k1 * h^15
err2 = k2 * h^7

integral error estimate = err1*(err1/err2)^2 = k1*(k1/k2)^2 * h^(15+8*2) = k * h^31

If (k, k1, k2) are similar in size, constant term matched too.

Exponents cannot be picked in random; it had to produce O(h^31) on the right.
If we use 5 or 7 nodes, instead of square, it would need some non-integer exponents.

Okay, I must be missing something. Why does it have to produce O(h^31) on the right? What if the algorithm turned out to be O(h^33) or some other order?

How would using a different number of nodes produce non-integer exponents?

Even if it did produce non-integer exponents, what would be wrong with that?

Sorry if I'm missing something obvious.
Find all posts by this user
Quote this message in a reply
08-01-2022, 07:01 PM
Post: #27
RE: Numerical integration methods
The exponent of the error estimate must match the error majoration of the best quadrature.
Find all posts by this user
Quote this message in a reply
08-01-2022, 07:24 PM
Post: #28
RE: Numerical integration methods
(08-01-2022 07:01 PM)parisse Wrote:  The exponent of the error estimate must match the error majoration of the best quadrature.

Okay, that's the piece I was missing. Everything else falls into place now.
Thank you.
Find all posts by this user
Quote this message in a reply
08-04-2022, 12:54 PM
Post: #29
RE: Numerical integration methods
(07-19-2022 05:16 PM)parisse Wrote:  int uses an adaptive method, as described here in French (Hairer: https://www.unige.ch/~hairer/poly/chap1.pdf).
Source code is available from giac source code
https://github.com/geogebra/giac/blob/ma...pp/intg.cc, function
Code:
  bool tegral(const gen & f,const gen & x,const gen & a_,const gen &b_,const gen & eps,int nmax,gen & value,bool exactcheck,GIAC_CONTEXT){...}

Bernard, out of interest, is a version of the Risch Algorthm still used for symbolic integration in XCAS/HP Prime?
Find all posts by this user
Quote this message in a reply
08-05-2022, 09:45 AM
Post: #30
RE: Numerical integration methods
Yes, the rational version is implemented.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: