Numerical integration methods
07-25-2022, 10:36 PM
Post: #21
 lrdheat Senior Member Posts: 856 Joined: Feb 2014
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!
07-26-2022, 05:18 PM
Post: #22
 Nigel (UK) Senior Member Posts: 518 Joined: Dec 2013
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)
07-31-2022, 04:59 PM
Post: #23
 parisse Senior Member Posts: 1,311 Joined: Dec 2013
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
08-01-2022, 03:43 AM
Post: #24
 Wes Loewer Senior Member Posts: 453 Joined: Jan 2014
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.
08-01-2022, 03:01 PM
Post: #25
 Albert Chan Senior Member Posts: 2,565 Joined: Jul 2018
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.
08-01-2022, 03:46 PM
Post: #26
 Wes Loewer Senior Member Posts: 453 Joined: Jan 2014
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.
08-01-2022, 07:01 PM
Post: #27
 parisse Senior Member Posts: 1,311 Joined: Dec 2013
RE: Numerical integration methods
The exponent of the error estimate must match the error majoration of the best quadrature.
08-01-2022, 07:24 PM
Post: #28
 Wes Loewer Senior Member Posts: 453 Joined: Jan 2014
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.
08-04-2022, 12:54 PM
Post: #29
 jonmoore Senior Member Posts: 309 Joined: Apr 2020
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?
08-05-2022, 09:45 AM
Post: #30
 parisse Senior Member Posts: 1,311 Joined: Dec 2013
RE: Numerical integration methods
Yes, the rational version is implemented.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)