Numerical integration methods

[Tonig00]

Hello

I would like to say that the prime is great.

Eventhough, sometimes there is some problem with a function. This is the case with the following integral (sorry, I do not know how to put a screenshot):
Int((x^2*((e^(1/x)^3)-1),x,1,100)

This, performed with a G2 with ROM 14603, takes around a minute giving 4.805 (in home mode).

I think it is a long time for this hardware.

Does anybody knows which numerical method the Prime uses?

Thanks very much

Toni[/quote]

[Albert Chan]

When I tried this, integral gives an warning "Adaptive method failure, will try with Romberg ..."
So, the trick is to "adapt" manually, by splitting the integral.

CAS> f(x) := x^2 * expm1(1/x^3)
CAS> int(f(x), x, 1., 20.), int(f(x), x, 20., 100.)

[3.19558488078, 1.60945857854]

∫ + ∫ = 4.80504345932

Udpate:

I increased HOME eps setting to 1e-9 (default was 1e-12)

CAS> int(f(x), x, 1., 100.)

4.80504345904

---

I had translated robve's double exponential quadrature code for HP Prime, here

CAS> quad(f, 1, 100)

[4.80504346031, 4.01306988706e−11]

quad(f,a,b) was defaulted with eps = 1e-9, thus slight error on last digit.

CAS> quad6(f, 1, 100, 1, 6, 1e-12)

[4.80504346032, 1.28096166025e−13]

[KeithB]

I don't know which method the int function uses, but you can force it to use romberg by using the romberg() command. (I suspect that int uses romberg when you ask for the definite integral).

[Albert Chan]

With substitution, HP Prime can solve integral symbolically.

CAS> f(x) := x^2 * expm1(1/x^3)
CAS> r := quote(int(f(x), x, a, b)) (x = y^(-1/3))

Touched up ⇒ \(\displaystyle {1\over3}\int _\frac{1}{a^3} ^\frac{1}{b^3} \frac{1-e^y}{y^2}\;dy\)

CAS> r := simpilfy(r)

1/3*(-a^3*e^(1/a^3) + b^3*e^(1/b^3) + Ei(1/a^3) - Ei(1/b^3) + a^3 - b^3)

CAS> r(a=1., b=100.)      → 4.80504345652

This is very bad numerical answer!
Let's replace exp with more accurate expm1 ... carefully!

CAS> g(x) := (x^3*expm1(1/x^3) - Ei(1/x^3))/3
CAS> g(100.) - g(1.)       → 4.80504346032

We write g as function, to prevent expm1(x) from expanding to exp(x)-1
If argument for g is numerical, expm1() is used, instead of exp()

[parisse]

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){...}

[lrdheat]

Just want to thank parisse for all the work that he contributed to the Prime. It looks as if there will not be much more attention given to the Prime by HP or Royal…you really did amazing work in the time that you had when HP had some interest in this, and I am glad to have this device.

[Albert Chan]

(07-18-2022 06:51 PM)Tonig00 Wrote:  Int((x^2*((e^(1/x)^3)-1),x,1,100)

This, performed with a G2 with ROM 14603, takes around a minute giving 4.805 (in home mode).

I think it is a long time for this hardware.

Does anybody knows which numerical method the Prime uses?

Turns out the issue is less related to integration routine.
It is more related to inability to evaluate integrand accurately.
Without this fuzziness, Gaussian Quadrature should not have trouble.

Rewrite expression using expm1() does not help. It get "simplified" away.
We might as well split up integrals, with both ∫ accurately calculated.

CAS> ∫(x^2*exp(1/x^3), x, 1., 100.), ∫(-x^2, x, 1., 100.)

[333337.805043, -333333.]

CAS> sum(Ans)

4.80504345335

log2(333333) ≈ 18.3 bits ≈ 5.5 decimal digits.
Decimal accuracy ≈ 14 - 5.5 = 8 to 9 digits

---

If expm1() work as advertised, apply x = exp(y) substitution is better.
This is because integrand looks very similar to 1/x = ln(x)'

∫(x^2*expm1(1/x^3),x,1.,100.) = ∫(exp(3y)*expm1(1/exp(3y)),y,0.,ln(100.))

lua> Q = require'quad'
lua> f = function(x) return x^2 * expm1(1/x^3) end
lua> g = function(x) x=exp(3*x); return x*expm1(1/x) end

It is difficult to fit 1/x well as polynomial.
For example, Q.gauss() use n=20 Gaussian Quadrature Weights and Abscissae

lua> Q.gauss(f, 1, 100)
4.778930417963698
lua> Q.gauss(g, 0, log(100))
4.805043460316733

lua> Q.gauss(f, 1, 10) + Q.gauss(f, 10, 100)
4.805043229893461
lua> Q.gauss(g, 0, log(10)) + Q.gauss(g, log(10), 2*log(10))
4.805043460319849

The transformation also helped tanh-sinh quadrature method
Note: Q.quad() have default eps = 1e-9

lua> Q.quad(f, 1, 100)
4.805043460309022         4.0108968271012166e-011     111
lua> Q.quad(g, 0, log(100))
4.8050434603198475       1.1409030351596591e-009     58

[Albert Chan]

(07-18-2022 09:26 PM)Albert Chan Wrote:  CAS> f(x) := x^2 * expm1(1/x^3)
CAS> r := quote(int(f(x), x, a, b)) (x = y^(-1/3))

Touched up ⇒ \(\displaystyle {1\over3}\int _\frac{1}{a^3} ^\frac{1}{b^3} \frac{1-e^y}{y^2}\;dy\)

(1-e^y)/(y*y) = -1/y - 1/2! - y/3! - y^2/4! - ...

Remove -1/y, integrand is a simple polynomial, easy to integrate.
Bonus, correction term is small; rough estimate is adequate.

1/3*∫((1-exp(y))/(y*y), y, 1/a^3, 1/b^3) = ln(b/a) - 1/3*∫((exp(y)-1-y)/(y*y), y, 1/a^3, 1/b^3)

NOTE: HP71B uses Romberg's method, with function u-transformed.

10 A=1 @ B=100 @ P=1E-6
20 DEF FNF(X) @ N=N+1 @ FNF=X^2*EXPM1(1/X^3) @ END DEF
30 DEF FNG(X) @ N=N+1 @ X=EXP(3*X) @ FNG=X*EXPM1(1/X) @ END DEF
40 DEF FNH(X) @ N=N+1 @ FNH=(EXPM1(X)-X)/(X*X) @ END DEF
50 N=0 @ DISP INTEGRAL(A,B,P,FNF(IVAR)),N
60 N=0 @ DISP INTEGRAL(LN(A),LN(B),P,FNG(IVAR)),N
70 N=0 @ DISP LN(B/A) - INTEGRAL(1/A^3,1/B^3,P,FNH(IVAR))/3,N

>RUN
4.80504346442      511
4.80504345776      127
4.80504346032      31

>P=1E-10
>RUN 50
4.80504346032      2047
4.80504346032      255
4.80504346032      63

[Albert Chan]

(07-21-2022 03:12 PM)Albert Chan Wrote:  Rewrite expression using expm1() does not help. It get "simplified" away.

It would be nice if user requested expm1(), expm1() is what you get.
As a patch, we can make our own version, that does not get simplify away.

(04-02-2014 05:18 PM)Thomas Klemm Wrote:  You could use \(e^x-1=2 \cdot sinh(\frac{x}{2}) \cdot e^{\frac{x}{2}}\) instead.

Code:

#cas
exp1(z) := 2*sinh(z/2)*exp(z/2);

log1(z) :=
BEGIN
  local y, z2;
  y := ln(1+z);
  z2 := exp1(y);
  RETURN y + (z-z2)/(z2+1); // taylor correction
END;
#end

log1(x) = log(1+x); exp1(x) = exp(x)-1
log1(x) is coded using result of exp1(), again to avoid it simplified away.

Bonus, this version work for complex numbers

Quote:Turns out the issue is less related to integration routine.
It is more related to inability to evaluate integrand accurately.
Without this fuzziness, Gaussian Quadrature should not have trouble.

Confirmed.

CAS> int(x*x*exp1(1/x^3), x, 1., 100.)                              → 4.80504346032
CAS> int(exp(3y)*exp1(exp(-3y)), y, 0., ln(100))                 → 4.80504346032
CAS> ln(100/1) - int((exp1(x)-x)/(x*x),x,1.,1/100.^3)/3      → 4.80504346032

[lrdheat]

Equation, in original form, integrates quickly to correct approximation on CASIO CLASSPAD 400. The TI Nspire ii (non cas) , if I am keying it in correctly comes up with an answer of >99. My 34S on a DM 42, if I am keying it in properly, produces an answer of 99+ as well.

[lrdheat]

TI Nspire CAS ii also produces 99+ answer if I am keying it in correctly…

[lrdheat]

Wow! My HP 35S came up with the correct answer on fix 4 in 45 seconds!

[Nigel (UK)]

(07-23-2022 08:12 PM)lrdheat Wrote:  TI Nspire CAS ii also produces 99+ answer if I am keying it in correctly…

I think the problem with the "99+" answers comes from the expression e^(1/x)^3. As written it's ambiguous, giving different answers depending on whether the powers are evaluated left-to-right (\({\rm e}^{3/x}\)) or right-to-left (\({\rm e}^{1/x^3}\)). The WP34S gives the correct answer when provided with the correct input, and I expect the TI machines will too.

Nigel (UK)

[lrdheat]

I am floored by this…The CASIO 9750 giii, with the equation entered precisely as in the CLASSPAD 400 produces the 99+ answer. It is in math in settings. Same in math-mix mode. Any monkeying around with parentheses is to no avail.

I will mess around with my WP 34S to see what I am doing wrong there. Glad to know that it can be done, I am probably overlooking something simple!

[Wes Loewer]

(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).

A huge thanks to Parisse for adding this adaptive method to the Prime. The Romberg method made sense back in the day when calculators had very limited memory, but with ample memory now available, adaptive methods are more appropriate. Certain integrands with cusps or asymptotes could slow Romberg to a crawl. The Prime usually makes short work of these.

Since many of us don't read French (myself included), allow me to summarize the article. (Some years ago with the help of Google Translate, I read through this article and I think I've got the gist of it. Corrections welcome.)

The Prime's adaptive method is a variation on Gaussian quadrature. Many are familiar with the more commonly used Gauss-Kronrod Method, so let me use that as a starting point for comparison. Calculators typically use a 15-node Gauss-Kronrod implementation that starts out with a 7-node Gaussian quadrature (exact for polynomials up to degree 13). Then, those 7 nodes are reused along with 8 additional Kronrod nodes using different weights for all 15 nodes to get a more accurate result (exact up to degree 22). The difference between these two calculations is used to estimate the error.

The Prime's method described in the Hairer document starts off with a 15-node Gaussian quadrature (exact up to degree 29). The integral is recalculated a 2nd time using 14 of the original 15 Gaussian nodes with different weights (exact up to degree 13). The integral is then recalculated a 3rd time using 6 of the original Gaussian nodes, again with different weights (exact up to degree 5). These three results are used to estimate the error.

In both methods, the function is evaluated a total of 15 times. If the estimated error is small enough, the more accurate calculation is used. If the error is too large, then the region is cut in half and each half is recursively calculated.

(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.)

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.]

[Nigel (UK)]

(07-24-2022 03:32 PM)lrdheat Wrote:  I am floored by this…The CASIO 9750 giii, with the equation entered precisely as in the CLASSPAD 400 produces the 99+ answer. It is in math in settings. Same in math-mix mode. Any monkeying around with parentheses is to no avail.

I will mess around with my WP 34S to see what I am doing wrong there. Glad to know that it can be done, I am probably overlooking something simple!

I don't have exactly these two calculators, but I do have a Casio 9860gii and a Casio Classpad 300+.

On the Casio 9860gii:

• In Math input mode, pressing 2^3^4 displays \(2^{3^4}\) and evaluates to \(2^{81}\).
• In Linear input mode, 2^3^4 displays 2^3^4 and evaluates to 4096.
  

On the Classpad 300+, pressing 2^3^4 displays 2^3^4 and evaluates to \(2^{81}\). I don't think that there is a math input mode on the 300+.

The program I'm using on my WP34S is LBL A x^2 x<>y 1/x 3 y^x e^x 1 - * RTN.

Nigel (UK)

[lrdheat]

Hi Nigel, I have not tried the WP 34S program, but graphing the function on the CLASSPAD 400, for x=2, f(x)=~.5326, same on HP 35S. Using x=2, and not using a program, just stepping through the equation manually using the HP 42S with it’s RPN (which I love), I cannot produce anything resembling.5326! I’ve come to the conclusion that I do not understand this function at all!

[Albert Chan]

(07-25-2022 03:27 AM)lrdheat Wrote:  Using x=2, and not using a program, just stepping through the equation manually using the HP 42S with it’s RPN (which I love), I cannot produce anything resembling.5326! I’ve come to the conclusion that I do not understand this function at all!

f(x) = x^2 * expm1(1/x^3), evaluated at x=2

2
ENTER X^2 *
LASTX X<>Y 1/X             // (1/x^3) (x^2)
E^X 1 - *

0.532593812267305267316028911247176

Or, more accurately, "E^X-1" for E^X 1 -

0.5325938122673052673160289112471756

[lrdheat]

Hi Nigel, Your program produces the desired result, and the term produces the .5326 result when using x=2. What is hanging me up is it appears as if you are multiplying the result of x^2 by (e^((1/x)^3)-1) instead of (e^((1/x)^3)_1) being multiplied by the 2, and then raising “x” to that product. The equation did not appear to be ambiguous to me, but I guess it is open to more than one interpretation!

[lrdheat]

(e^((1/x)^3)-1), not _1 my typing mistake.