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 Code: This is a very difficult integral to compute numerically to any decent accuracy ... 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. |
|||
05-21-2021, 07:32 PM
Post: #2
|
|||
|
|||
RE: HP71 IBOUND fooled
Code: / Inf 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 |
|||
05-21-2021, 09:38 PM
Post: #3
|
|||
|
|||
RE: HP71B IBOUND fooled
(05-21-2021 07:32 PM)Albert Chan Wrote: F(∞) - F(0) 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 |
|||
05-02-2022, 01:42 AM
(This post was last modified: 05-03-2022 01:07 PM by Albert Chan.)
Post: #4
|
|||
|
|||
RE: HP71B IBOUND fooled
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\) // t0 = 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 https://en.wikipedia.org/wiki/Q-function |
|||
05-02-2022, 02:57 PM
Post: #5
|
|||
|
|||
RE: HP71B IBOUND fooled
(05-02-2022 01:42 AM)Albert Chan Wrote: Unlike Kahan's transformation, getting u end-point involved solving cubics. 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 |
|||
05-03-2022, 07:09 PM
(This post was last modified: 05-10-2022 11:04 AM by Albert Chan.)
Post: #6
|
|||
|
|||
RE: HP71B IBOUND fooled
(05-02-2022 01:42 AM)Albert Chan Wrote: 50 Y=EXP(-X*X/4) @ Y=Y/2/COS(ASIN(Y)/3+PI/6) 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 |
|||
08-10-2022, 04:48 PM
Post: #7
|
|||
|
|||
RE: HP71B IBOUND fooled
(05-02-2022 02:57 PM)Albert Chan Wrote: This solve cubics, x^3 = a*x + b, for x 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 |
|||
08-10-2022, 06:01 PM
(This post was last modified: 08-11-2022 01:12 PM by Albert Chan.)
Post: #8
|
|||
|
|||
RE: HP71B IBOUND fooled
(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 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)