HP71B IBOUND fooled

[Albert Chan]

Here is an integral, from a very old (2006) thread: So your HP can INTEGRATE ...

Code:

       / Inf
I1 =   |     sin(x)*sin(x^2) .dx
       / 0

Solution from Valentin:

Code:

This is a very difficult integral to compute numerically to any decent accuracy ...
Fortunately, there's a nifty closed form for a generalized family of similar integrals, 
this being but one of the simplest cases, namely:

                  / A                             / A 
    I1 = (Cos(A)* |   Cos(x)/Sqrt(x).dx + Sin(A)* |   Sin(x)/Sqrt(x).dx )/2
                  / 0                             / 0  
where A=1/4 and both non-elementary integrals are particular cases of Fresnel functions.

I tried in emu71, and showed timing and IBOUND

>A=1/4 @ P=1E-12
>X1=COS(A) @ Y1=SIN(A)

>T=TIME @ X2=INTEG(0,A,P,COS(IX)/SQRT(IX))/2 @ TIME-T, IBOUND
40.58       -9.93754844068E-13
>T=TIME @ Y2=INTEG(0,A,P,SIN(IX)/SQRT(IX))/2 @ TIME-T, IBOUND
.11           8.2963337046E-14
>X2, Y2, X1*X2+Y1*Y2
.496880004382       4.14810242686E-2       .491695777984

X2 required much time to calculate. But, IBOUND numbers seems good ...

However, if we rewrite integral, letting x=t^2, dx=2t dt, we get different numbers.
(we might as well do Y2 too, to compare effect of x=t^2 substitution)

>T=TIME @ X2=INTEG(0,.5,P,COS(IX*IX)) @ TIME-T, IBOUND
.11           4.96884313442E-13
>T=TIME @ Y2=INTEG(0,.5,P,SIN(IX*IX)) @ TIME-T, IBOUND
.17           4.14807134341E-14
>X2, Y2, X1*X2+Y1*Y2
.496884029215       4.14810242685E-2       .491699677694

We get exact result (all 12 digits correct !), using much less time.

---

X2 (with sqrt denominator), after u-transformation, end-point is not zero !
INTEGRAL assumed zero endpoints (both side), thus bad results.

\(\displaystyle \int_0^{1\over4} {\cos(x)\over2\sqrt{x}} dx
= \int_0^1 {\cos(t/4)\over4\sqrt{t}} dt
= \int_0^1 {6u(1-u)\cos(u^2(3-2u)/4)\over4\sqrt{u^2(3-2u)}} du
= {\sqrt{3}\over2}\int_0^1 \left(1-{2u\over3}-{u^2\over6}\;-\;...\right) du
\)

However, u-transformed plot looks like a straight line !
IBOUND numbers were fooled, suggesting excellent estimate, which it isn't.

[Albert Chan]

Code:

       / Inf
I1 =   |     sin(x)*sin(x^2) .dx
       / 0

We can solve above integral, using erf():

Let cis(x) = exp(i*x)

∫(cis(x^2))
= ∫(exp(-(w*x)^2)), where w = cis(-pi/4)
= sqrt(pi)/2 * erf(x*w)/w

∫(cis(x^2±x))
= ∫(cis((x±1/2)^2 - 1/4))
= cis(-1/4) * ∫(cis((x±1/2)^2))
= cis(-1/4) * sqrt(pi)/2 * erf((x±1/2)*w)/w

F(x) = ∫(sin(x)*sin(x^2)) = ∫(sinh(i*x)/i * sinh(i*x^2)/i)

= -1/4 * ∫((cis(x)-cis(-x)) * (cis(x^2)-cis(-x^2)))
= -1/4 * ∫((cis(x+x^2)+cis(-x-x^2)) - (cis(-x+x^2)+cis(x-x^2)))
= -1/2 * re(∫(cis(x^2+x) - cis(x^2-x)))
= -1/2 * re(cis(-1/4) * sqrt(pi)/2 * (erf((x+1/2)*w) - erf((x-1/2)*w))/w)

erf(∞) = 1 -> F(∞) = 0

erf() is odd -> F(0) = - re(cis(-1/4) * sqrt(pi)/2 * erf(w/2)/w)

Power Series Expansion of the Error Function
sqrt(pi/2)*erf(z) = z - z^3/(3*1!) + z^5/(5*2!) - z^7/(7*3!) + ...

sqrt(pi)/2 * erf(w/2)/w
= 1/2                - (-i)/(3*1!*2^3)     + (-1)/(5*2!*2^5)    - (+i)/(7*3!*2^7)
+ 1/(9*4!*2^9) - (-i)/(11*5!*2^11) + (-1)/(13*6!*2^13) - (+i)/(15*7!*2^15) + ...
= (1/2 - 1/320 + 1/110592 - 1/76677120 + ...)
+ (1/24 - 1/5376 + 1/2703360 - 1/2477260800 + ...) * i
≈ 0.496884029215 + 0.0414810242685*i

F(∞) - F(0)
= re(cis(-1/4) * sqrt(pi)/2 * erf(w/2)/w)
≈ cos(1/4) * 0.496884029215 + sin(1/4) * 0.0414810242685
≈ 0.491699677694

[Albert Chan]

(05-21-2021 07:32 PM)Albert Chan Wrote:  F(∞) - F(0)
= re(cis(-1/4) * sqrt(pi)/2 * erf(w/2)/w)
≈ cos(1/4) * 0.496884029215 + sin(1/4) * 0.0414810242685
≈ 0.491699677694

Just found a more elegant way, without doing even cos(1/4), sin(1/4).

Abramowitz and Stegun, eqn 7.1.6: \(\quad\displaystyle erf(z) =
{2\over\sqrt{\pi}} e^{-z^2} \sum_{n=0}^∞ {2^n \over 1·3· \cdots (2n+1)} z^{2n+1} \)

exp(-z^2) = exp(-w^2/4) = exp(i/4) = cis(1/4)
This got completely cancelled with term cis(-1/4)

This left only summation terms ... and only the real terms

F(∞) - F(0)
= 1/2 - 1/(3*5*2^3) + 1/(3*5*7*9*2^5) - 1/(3*5*7*9*11*13*2^7) + ...
= 1/2 * (1 - 1/(3*5*4) * (1 - 1/(7*9*4) * (1 - 1/(11*13*4) * (1 - ...))))
≈ 0.491699677694

[Albert Chan]

From Kahan's https://people.eecs.berkeley.edu/~wkahan/MathSand.pdf, page 23

\(\displaystyle Q(x) = \frac{1}{\sqrt{2\pi}} \int_x^\infty \exp\left(-\frac{u^2}{2} \right)\;du \)

To make integral limits finite, we let t = exp(-u*u/2) ⇒ dt = (t) (-u) du

\(\displaystyle Q(x) = \frac{1}{\sqrt{2\pi}}
\int_0^{t_0} \frac{dt}{\sqrt{-2 \ln(t)}} \qquad\) // t₀ = exp(-x²/2) ≤ 1

For HP71B, this integral have problems: t=0, infinite slope; t=1, infinite value.
INTEGRAL internal u-transform is enough to make end-points finite, but not zero.

Kahan's method "B": Let t = sin(s²) ⇒ dt/ds = cos(s²) (2s)
For both t=0 (s=0) or t=1 (s=√(pi/2)), dt/ds = 0

We could also work in reverse, with dt/du = (1-u) (6u) ⇒ t = u² (3-2u)
Unlike Kahan's transformation, getting u end-point involved solving cubics.
Luckily, with 0 ≤ x ≤ 1, only 1 u solution satisfied.

\(\displaystyle \large \int_0^x f(t)\;dt \qquad\) // 0 ≤ x ≤ 1

\(\displaystyle = \int_0^{s_0} f(\; \sin(s^2)\; )\;\cos(s^2)\;(2s)\;ds\qquad\quad\) // s0 = sqrt(asin(x))
\(\displaystyle = \int_0^{u_0} f(3u^2-2u^3)\;(1-u)\;(6u)\;du\qquad\,\) // u0 = (sqrt(x)/2) / cos(asin(sqrt(x))/3 + pi/6))

We could also scale integral to limits 0 to 1, then do u-transform. (u0 = 1)
But, this left scaling factor inside integrand, not as efficient.
Getting u0 involve only 1 time cost.

10 INPUT "x,eps = ? ";X,P
20 Y=SQR(ASIN(EXP(-X*X/2)))
30 DEF FNS(T) @ S=SIN(T*T) @ FNS=T*SQR((S+1)*(S-1)/LN(S)) @ END DEF
40 T=TIME @ DISP INTEGRAL(0,Y,P,FNS(IVAR))/SQR(PI),TIME-T

50 Y=EXP(-X*X/4) @ Y=Y/2/COS(ASIN(Y)/3+PI/6)
60 DEF FNU(U)=U*(1-U)/SQR(-LN(U*U*(3-2*U)))
70 T=TIME @ DISP INTEGRAL(0,Y,P,FNU(IVAR))*3/SQR(PI),TIME-T
80 GOTO 10

Code:

>RUN
x,eps = ? 0, 1E-3
 .500003367483        .44
 .500002347538        .37
x,eps = ? 0, 1E-6
 .500000000524        1.52
 .500000057737        .7
x,eps = ? 0, 1E-9
 .499999999995        5.79
 .500000004166        1.88

x,eps = ? 10, 1E-3
 7.61991738243E-24    .21
 7.61991738248E-24    .23
x,eps = ? 10, 1E-6
 7.61985318691E-24    .93
 7.61985318696E-24    .67
x,eps = ? 10, 1E-9
 7.61985302414E-24    5.83
 7.61985302419E-24    4.87

x,eps = ? 20, 1E-3
 2.75363048487E-89    .19
 2.75363048486E-89    .21
x,eps = ? 20, 1E-6
 2.75362414276E-89    .62
 2.75362414274E-89    .7
x,eps = ? 20, 1E-9
 2.75362411868E-89    2.83
 2.75362411866E-89    2.25

X,EPS = ? 30, 1E-3
 4.90671903956E-198   .22
 4.90671903951E-198   .24
X,EPS = ? 30, 1E-6
 4.90671417978E-198   .52
 4.90671417974E-198   .42
X,EPS = ? 30, 1E-9
 4.90671392717E-198   2.85
 4.90671392713E-198   2.4

https://en.wikipedia.org/wiki/Q-function

[Albert Chan]

(05-02-2022 01:42 AM)Albert Chan Wrote:  Unlike Kahan's transformation, getting u end-point involved solving cubics.
Luckily, with 0 ≤ x ≤ 1, only 1 u solution satisfied.

This solve cubics, x^3 = a*x + b, for x
Formula is general for (a,b), a ≠ 0, even complex numbers.
Note: constant k is to transform cubic to form, sin(3θ) = 3*sin(θ) - 4*sin(θ)^3

XCas> cubic(a,b) := {local k:=sqrt(4*a/3); k .* sin(([0,2*pi,-2*pi].-asin(4b/k^3))./3)}

3u² - 2u³ = x
3/u - 2 = x/u³
(1/u)³ = (3/x)(1/u) + (-2/x)

XCas> simplify(1 ./ cubic(3/x, -2/x))

\(\displaystyle \left[
\frac{\sqrt{x}}{2 \sin\left(\frac{\mathrm{asin}\left(\sqrt{x}\right)}{3}\right)},
\frac{\sqrt{x}}{2 \sin\left(\frac{\mathrm{asin}\left(\sqrt{x}\right)+2\cdot \pi}{3}\right)},
\frac{\sqrt{x}}{2 \sin\left(\frac{\mathrm{asin}\left(\sqrt{x}\right)-2\cdot \pi}{3}\right)}
\right]\)

Consider only denominators, with asin(√x)/3 = 0 .. pi/6
Denominator = 2*sin([(0 .. pi/6), (2/3*pi .. 5/6*pi), (-2/3*pi .. -pi/2)]) = [ (0 .. 1), (1..√3), (-2..-√3) ]

Only middle root gives correct u range:

3u² - 2u³ = x ≥ 0
u² = x/(3-2u) ≤ x/(3-2*max(u)) = x/1
0 ≤ u ≤ √x

[Albert Chan]

(05-02-2022 01:42 AM)Albert Chan Wrote:  50 Y=EXP(-X*X/4) @ Y=Y/2/COS(ASIN(Y)/3+PI/6)
60 DEF FNU(U)=U*(1-U)/SQR(-LN(U*U*(3-2*U)))
70 T=TIME @ DISP INTEGRAL(0,Y,P,FNU(IVAR))*3/SQR(PI),TIME-T

What if x is big ? (upper integral limit, u0 is tiny. Or, 1-u ≈ 1)

\(\displaystyle Q(x)
= {3 \over \sqrt{\pi}} \int_0^{u_0} \frac{u\;(1-u)\;du}{\sqrt{-\log(u^2\;(3-2u))}}
≈ {3 \over \sqrt{2\;\pi}} \int_0^{u_0} \frac{u\;du}{\sqrt{-\log(u)}}
= 3 \; Q\left(\sqrt{-4\;\log(u_0)}\;\right)
\)

XCas> y := exp(-x*x/4)
XCas> u0 := y/2 / cos(asin(y)/3 + pi/6)
XCas> series(2*sqrt(-log(u0)), x=inf, polynom)      → \(\displaystyle x + \frac{\ln\left(3\right)}{x}\)

For big x, we have Q(x) / Q(x + ln(3)/x) ≈ 3
Perhaps we can replace constant 3 by other multiplier M ?

\(\displaystyle M =
\frac{Q(x)}{Q \left(x+\frac{k}{x} \right)}
= \frac{erfc \left({x \over \sqrt2}\right)/2} {erfc \left({x \over \sqrt2}
+ {k \over \sqrt2 \; x} \right)/2}
= e^k \left[ 1
+ \frac{k\;(k+2)}{2 x^2}
+ \frac{k\;(k^3+4k^2-16)}{8x^4} \;+\; ... \right]\)

Example: find y such that Q(y) = Q(x = 5) / 2

lua> x, m = 5, log(2)
lua> k = m
lua> x + k/x -- rough y
5.138629436111989
lua> k = m - log1p(k*(k+2)/(2*x^2))
lua> k = m - log1p(k*(k+2)/(2*x^2))
lua> x + k/x
5.131772468987268
lua> k = m - log1p(k*(k+2)/(2*x^2) + k*(k^3+4*k^2-16)/(8*x^4))
lua> k = m - log1p(k*(k+2)/(2*x^2) + k*(k^3+4*k^2-16)/(8*x^4))
lua> x + k/x
5.132077991427134

lua> -icdf(cdf(-x)/2) -- actual y
5.132018332044298

Comment: For k correction, we do not need accurate log1p()
Even simple log1p(ε) ≈ ε may suffice.

Or, we can revert series for k

lua> t1 = - m*(m+2)/(2*x^2)
lua> t2 = m*(m^2+4*m+6)/(2*x^4)
lua> t3 = - m*(15*m^3+92*m^2+252*m+360)/(24*x^6)
lua> k = t3 + t2 + t1 + k
lua> x + k/x
5.131972797544282

To improve y estimate further, apply Aitken Delta^2 process

lua> k = k - t3*t3 / (t3 - t2)
lua> x + k/x
5.132010307411825

---

9/22/2024, Simpler proof, based from Q continued fraction

\(\displaystyle Q(x ≥ 0) =
\frac{\text{pdf (x)}}{x +}\;
\frac{1}{x+}\; \frac{2}{x+}\; \frac{3}{x+}\;...
\)

If x is huge, we can estimate with just first CF "term"

\(\displaystyle M =
\frac{Q(x)}{Q \left(x+\frac{k}{x} \right)}
≈\left( \frac
{\frac{1}{\sqrt{2\;\pi}} \exp(-\frac{x^2}{2}) ÷ x}
{\frac{1}{\sqrt{2\;\pi}} \exp(-\frac{(x+\frac{k}{x})^2}{2}) ÷ (x+\frac{k}{x})}
\right)
= \exp\left(k + \frac{k^2}{2\;x^2}\right) \left(1 + \frac{k}{x^2}\right)
≈ e^k
\)

[Albert Chan]

(05-02-2022 02:57 PM)Albert Chan Wrote:  This solve cubics, x^3 = a*x + b, for x
Formula is general for (a,b), a ≠ 0, even complex numbers.
Note: constant k is to transform cubic to form, sin(3θ) = 3*sin(θ) - 4*sin(θ)^3

XCas> cubic(a,b) := {local k:=sqrt(4*a/3); k .* sin(([0,2*pi,-2*pi].-asin(4b/k^3))./3)}

For cubic(a,b), (sin,asin) → (cos,acos), we *still* get solutions of cubic, albeit in different order.

XCAS> cubic2(a,b) := {local k:=sqrt(4*a/3); k .* cos(([0,2*pi,-2*pi].-acos(4b/k^3))./3)}

If (a,b) are real, complex roots, if any, always comes in conjugate pairs.
With above symmetry, if asin(), acos() are real, this implied at least 2 real roots --> all cubic roots real.

(cos, acos) can easily mapped to (cosh, acosh)

cosh(z) = cos(±i*z)            // both cosh and cos are even function.
acosh(z) = acos(z) * ±i      // pick sign to produce non-negative real part

Example, for value of Super Golden Ratio

x^3 = x^2 + 1
(x-1/3)^3 = (x-1/3) * (1/3) + (29/27)      // depressed cubic

(x-1/3) = 2/3 * cos((2*n*pi - acos(29/2))/3)

29/2 > 1, cubic had only 1 real root, for n=0:

x = 2/3*cosh(acosh(29/2)/3) + 1/3 ≈ 1.465571232

[Albert Chan]

(08-10-2022 04:48 PM)Albert Chan Wrote:  For cubic(a,b), (sin,asin) → (cos,acos), we *still* get solutions of cubic, albeit in different order.

Let θ = asin(4b/sqrt(4*a/3)^3)

(0 - θ)/3 = pi/2 - (2*pi - (pi/2-θ))/3
(2*pi - θ)/3 = pi/2 - (0 - (pi/2-θ))/3

Swapping sin/asin to cos/acos, or vice versa, first 2 solutions order also swapped.

Again, Super Golden ratio example, for all roots.

XCAS> cubic(1/3, 29/27) .+ 1/3.

[-0.232786+0.792552*i, 1.46557+5.60374e-17*i, -0.232786-0.792552*i]

XCAS> cubic2(1/3, 29/27) .+ 1/3.

[1.46557, -0.232786+0.792552*i, -0.232786-0.792552*i]

---

Explanation of algebraic way to solve cubic roots
https://www.hpmuseum.org/forum/thread-10...#pid150296

Kahan's method, for accurate numerical roots of cubic
https://math.stackexchange.com/a/3374897