Post Reply 
HP71B IBOUND fooled
05-21-2021, 07:17 PM (This post was last modified: 05-21-2021 10:46 PM by Albert Chan.)
Post: #1
HP71B IBOUND fooled
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.
Find all posts by this user
Quote this message in a reply
05-21-2021, 07:32 PM
Post: #2
RE: HP71 IBOUND fooled
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
Find all posts by this user
Quote this message in a reply
05-21-2021, 09:38 PM
Post: #3
RE: HP71B IBOUND fooled
(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 Smile

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
Find all posts by this user
Quote this message in a reply
Post Reply 




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