arcsinc( 1-y ), for small y
07-05-2018, 11:43 PM
Post: #1
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
arcsinc( 1-y ), for small y
Come across a fun civil engineering problem from HP archive [1]
"Five minute challenge" essentially evaluate arcsinc( 1-y )

sin(x) / x = ( 1-y ) -- y = 1/5281 = rail expansion effect
x = arcsinc( 1-y ) -- then, rail warp e = 2640 tan(x/2)

From Gerson W. Barbosa (msg #21), sinc(x) ~ cos(x / sqrt(3))

x = arcsinc( 1-y ) ~ sqrt(3) arccos( 1-y )
= 0.03370733324 -- using my Casio fx-115ms

Estimate is good, exact x = 0.033707758809879429 ...

I wanted a simpler, and more accurate arcsinc( 1-y )

arccos( 1-y ) taylor series = sqrt(2y) (1 + y/12 + ...) [2]

x = arcsinc( 1-y ) ~ sqrt(6y) = 0.03370680134

For better accuracy, apply Newton's method, x0 = sqrt(6y):

f(x) = sinc(x) - (1-y)
= y - x^2/6 + x^4/120 - x^6/5040 + ...
f'(x) = -x/3 + x^3/30 - (3/70) x^5 + ...

f(x0) = y - y + 0.3 y^2 - (3/70) y^3 + ... ~ 0.3 y^2
f'(x0) = x0 (-1/3 + y/5 - (54/35) y^2 + ...) ~ -x0 / 3

x = x0 - f(x0) / f'(x0)
~ x0 + (0.9 y^2 / x0) * (x0/x0)
= x0 (1 + 0.9 y^2 / (6y))
= sqrt(6y) (1 + 0.15 y)

x = arcsinc( 1-y ) ~ sqrt(6y) (1 + 0.15 y) = 0.03370775874

with x this good, f(x) = 0.0 on my Casio :-)

[2] https://math.stackexchange.com/questions...-two-terms
07-06-2018, 03:22 PM (This post was last modified: 10-01-2019 10:54 PM by Albert Chan.)
Post: #2
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
This new formula is slightly more accurate, and just as simple:

For small y: arcsinc( 1-y ) ~ √(6y) / (1 - 0.15 y)

arcsinc(1 - 1/5281) ~ √(6/5281) / (1 - 0.15/5281) = 0.033707 75877

Example:
y = 1 - sin(x) / x = 0.001665833532

Old formula, x ~ √(6y) * (1 + 0.15 y) = 0.099999 98408
New formula, x ~ √(6y) / (1 - 0.15 y) = 0.099999 99033

Update: newer formula from post #11

x ~$$\large \sqrt{6y ÷ \sqrt{1 + 0.6 y}}$$ = 0.099999 99997 // almost round-trip!
07-06-2018, 09:07 PM
Post: #3
 Thomas Klemm Senior Member Posts: 1,804 Joined: Dec 2013
RE: arcsinc( 1-y ), for small y
(07-05-2018 11:43 PM)Albert Chan Wrote:  Come across a fun civil engineering problem from HP archive [1]

Thank you for reminding me of this thread.
Good memories.

Quote:"Five minute challenge" essentially evaluate arcsinc( 1-y )

And this was supposed to be a 5-minute challenge...

Well, probably not.

Cheers
Thomas
07-06-2018, 11:17 PM (This post was last modified: 07-10-2018 09:38 PM by Albert Chan.)
Post: #4
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
To show new formula is *really* more accurate: (below from Mathematica)

arcsinc( 1-y ) / sqrt(6y)
= 1 + 3/20 y + 321/5600 y^2 + 3197/112000 y^3 + 445617/27596800 y^4 + ...
~ 1 + 0.15 y + 0.0573214 y^2 + 0.0285446 y^3 + 0.0161474 y^4 + ...
~ 1 + 0.15 y

sqrt(6y) / arcsinc( 1-y )
= 1 - 3/20 y - 39/1120 y^2 - 1649/112000 y^3 - 5285631/689920000 y^4 - ...
~ 1 - 0.15 y - 0.0348214 y^2 - 0.0147232 y^3 - 0.00766122 y^4 - ...
~ 1 - 0.15 y

Looking at the size of dropped terms, new formula is more accurate:

arcsinc( 1-y ) ~ sqrt(6y) / (1 - k y), where k = 0.15

For more accurate arcsinc, correct for k using above divide trick.

k = 0.15 / (1 - 13/56 y / (1 - 0.19068 y / (1 - 0.21627 y)))
07-07-2018, 06:11 AM (This post was last modified: 08-18-2019 01:14 PM by Albert Chan.)
Post: #5
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
I re-read the Five Minute challenge posts (highly recommended)
I found a simple formula that is slightly better than e = sqrt(6 m d) / 4

Let z = gain in length of of hypotenuse
Using post #30 3/4 rule, z = d/2 * 3/4 = 3/8 d

e^2 = (m/2 + z)^2 - (m/2)^2 = m z + z^2

e = sqrt((m + z) z)

With m=5280, d=1: e = sqrt((5280 + 3/8)(3/8)) = 44.49877105 ~ 44.50 ft

3/4 rule is the upper bound for tiny angle, so z (thus e) is over-estimated.

Let's try with bigger sector angle, say, d = 100 ft:

Old formula, e = sqrt(6 * 5280 * 100) / 4 = 444.9719092 ~ 444.97 ft
New formula, e = sqrt((5280 + 300/8)(300/8)) = 446.5492694 ~ 446.55 ft

We have a nice bound, correct e is between 444.97 ft to 446.55 ft

For correct value of e, with big angle, do arcsinc with corrected k:

y = 1 - 5280/5380 = 100/5380 = 0.01858736
k = 0.15 / (1 - 13/56 y / (1 - 0.19068 y / (1 - 0.21627 y))) = 0.1506523749
x = sqrt(6y) / (1 - k y) = 0.3348901066 (about 19 degree)

e = (m/2) tan( x/2 ) = 446.2332298 ~ 446.23 ft (correct)
07-07-2018, 06:04 PM (This post was last modified: 08-18-2019 01:17 PM by Albert Chan.)
Post: #6
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
Did plots of how good e = sqrt((m + z)z), vs e = sqrt(6 m d) / 4
The new formula is 4 times more accurate !

Because the ratio is so consistent around 4, this is a fast way to estimate e:
Previous example, d = 100 ft:

e ~ (446.5492654 * 4 + 444.9719092) / 5 = 446.2337942 ft

With 6.5X warp or less (d < 5.5 *miles*), accuracy can be improved:
y = 100/5380
w = 4 - 27/56 y = 3.991038237
e ~ (446.5492654 w + 444.9719092) / (w+1) = 446.2332309 ft

This is not far from exact e of 446.233229831394 ...

We could even reverse the process, use good e to estimate arcsinc !
To simplify, let m = 1, so e = 1/2 tan( x/2 )

Example: calculate arcsinc(0.99):

sinc(x) = 0.99 = 1/(1+ d)
d = 1/99
z = 3/8 d = 1/264

e(min) = √(z) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ = 0.06154574549
e(max) = √(z + z²) = 0.06166219923

w ~ 4 - 27/56 y = 3.995178571
e ~ (w e(max) + e(min)) / (w + 1) = 0.061638886
x ~ 2 arctan( 2 e ) = 0.2453178089

For comparison, asinc(0.99) = 0.2453178088 54025 ...
07-08-2018, 03:28 AM (This post was last modified: 07-13-2018 09:46 PM by Albert Chan.)
Post: #7
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
Redoing the Five minutes challenge with weighted e:
Note: sqrt(6 m d) / 4 = sqrt(6 m (8/3 z) / 16) = sqrt(m z)

m = 5280
z = 3/8 d = 3/8

e(min) = sqrt(m z) = 44.49719092 ~ 44.50
e(max) = sqrt((m + z) z) = 44.49877105 ~ 44.50

The nice error bound *guaranteed* 4 digits accuracy

e ~ (4 e(max) + e(min)) / 5 = 44.49845502

4 weighted e reach 10 digits accuracy !
w weighted e do even better, see next post ...
07-09-2018, 01:12 AM (This post was last modified: 07-13-2018 08:03 PM by Albert Chan.)
Post: #8
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
To show weighted e method work, I use Mathematica:
To simplify, assumed m = 1, so e = 1/2 tan(x/2)

sinc[x_] := Sin[x] / x;
t = Series[1 - sinc[Sqrt[x]], {x, 0, 5}];
t = InverseSeries[t] /. x -> y;
exact5 = 1/2 Tan[Sqrt[t] / 2];

w = 4 - 27/56 y;
weighted = (w * Sqrt[z + z z] + Sqrt[z]) / (w+1);
weighted = weighted /. z -> 3/8d /. d -> y/(1-y);

ratio = Series[exact5 / weighted, {y, 0, 5}] // Simplify // Normal
ratio // N

>> 1 - 9459/25088000 y^3 + 2651751/309084160000 y^4
>> 1. - 3.77033e-4 y^3 + 8.57938e-6 y^4

weighted e estimate is simple, and very good. :-)

To check if above is valid, use ratio to improve Five Minute Challenge:

First, calculate exact e, with 30 significant digits:

FindRoot[sinc[x] == 5280/5281, {x, 0.1}, WorkingPrecision -> 50]
>> {x -> 0.033707758809879429348218498305111 ... }

e = N[5280/2 Tan[x / 2] /. %, 30] (* Reference: all good digits *)
>> 44.4984550191007992545541600167

(* Normally no need for scaling, but above assumed m=1 *)

e = N[5280 weighted /. y-> 1/5281, 30] (* 14 digits accurate *)
>> 44.4984550191009131676665673849

e = N[5280 weighted ratio /. y-> 1/5281, 30] (* 22 digits accurate *)
>> 44.4984550191007992545475544533
07-12-2018, 05:26 PM (This post was last modified: 07-16-2018 01:54 PM by Albert Chan.)
Post: #9
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
On the railroad warp problem, same formulas still work for pi/2 < x < pi:
In other words, circle center and warp ok on the *same* side:

Prove sin(x) / x = m / (m + d):

m/2 = R sin(pi - x) -- chord length
(m + d) / 2 = R x -- arc length

sin(pi - x)/x = sin(x)/x = m / (m + d) -- QED

Prove e = (m/2) tan( x/2 ):

(m/2) / e = tan((pi - x) / 2) -- triangle: chord + e (contained circle center)

e = (m/2) / tan(pi/2 - x/2) = (m/2) tan( x/2 ) -- QED

My weighted e is not designed for big angle, but still work ok
For x < pi/2 (thus, warp < pi/2), weighted e has max error of 0.001%

For HUGE warp, weighted e formula run into problem:
Example: let m=1, 6.5X warp (d=5.5), calculate e:

sin(x) / x = 1 / (1+d) = 0.1538461538
x = 2.71131291 -- my Casio solver, x ~ 155 degrees
R = (1+d) / (2x) = 1.198681269
e = 1/2 tan(x/2) = 2.288101656 -- almost diameter of circle !

Assumed we don't know x, use 4 weighted e method:
z = 3/8 d = 2.0625
e(min) = sqrt(z) = 1.436140662
e(max) = sqrt(z + z z) = 2.513246158
e ~ (4 e(max) + e(min)) / 5 = 2.297825059 -- error = +0.425%

With warp that high, "improved" weight does not help:
y = 1 - 0.1538461538 = 0.8461538462
w = 4 - 27/56 y = 3.5920322967
e ~ (w e(max) + e(min)) / (w+1) = 2.27868654 -- error = -0.411%

For HUGE warp, weighted e method is not suited, below fit is better:
For pi/2 < x < pi, fitted e (below) has max error of 0.013%

e (warp 2.4X or more)
~ (((-0.690948 y + 1.95561) y - 2.47546) y + 1.52907) d

Redo above case (6.5X warp):

e ~ 0.4160261229 * 5.5 = 2.288143676 -- error = +0.002%
08-20-2018, 02:20 PM (This post was last modified: 08-21-2018 10:45 PM by Albert Chan.)
Post: #10
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
Forman Acton's book Numerical Method that work (page 68) have a neat solution to the railroad problem:

sinc(x) = 5280/5281 = 1/(1+e), where e=1/5280

--> 1 - x^2/3! + x^4/5! - x^6/7! + ... = 1 - e/(1+e)

Subtract 1 from both side, factor out x^2, we get:

x^2 = 6e / (1+e) / (1 - x^2/20 + ...)

Last denominator is about 1, we get x^2 ~ 6e / (1 + e) = 6/5281
We get x = sqrt(6/5281) = 0.0337068013

Iterate once more:

x^2 = 6e / (1 + e) / (1 - (6e / (1+e)) / 20) = 6e / (1 + 0.7e) = 6/5280.7
We get x = sqrt(6/5280.7) = 0.0337077588

With this x, warp = (5280/2) tan(x/2) = 44.498455 ft (matched rounded exact value)

The book use this example to illustrate two points:
1. Remove large masking quantities, (in this case, the 1.0 term)
2. Do not solve quadratic equation if the quadratic term carry little weight.

Amazing book ...
08-20-2018, 03:23 PM (This post was last modified: 06-18-2020 01:10 PM by Albert Chan.)
Post: #11
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
Rearange Acton fromula in terms of arcsinc(k = 1-y):

x = arcsinc(k) ~ sqrt(6y / (1 - 0.3 y))

this is much better than my estimate of sqrt(6y) / (1 - 0.15y)
For small y, errors cut down by almost 33% !

To push the idea a bit further, below cut down error (small y) by 97% !

x = arcsinc(k) ~ sqrt(6y / sqrt(1 - 0.6 y))

For y > 0.5, asymptotic formula is better: x ~ Pi / (1 + k + 1.244 k^3)

For 0<y<1, above combined setup kept relative error to 0.1% or less.

For more accuracy, use either form of Newton's method:

f = sin(x) - kx ﻿ ﻿ ﻿ ﻿ ﻿ ﻿→ x = x - (sin(x) - k*x) / (cos(x) - k)
f = sin(x)/x - k ﻿ ﻿ ﻿ ﻿ → x = x + x * (z - k) / (z - cos(x)), where z = sin(x)/x

Update: More accurate Asinc (small y), but more complicated formula:

Update2: instead of asymptotic formula, correction also work:
For y >= 0.4, asinc(1-y) ~ sqrt(6y / sqrt(1 - 0.6 y)) * (1.001 + y^6/53)
08-25-2018, 03:51 PM (This post was last modified: 10-02-2019 12:01 PM by Albert Chan.)
Post: #12
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
The book Numerical Method that work showed Acton's favorite way of doing interpolation (no difference table)
A variation of Aitken scheme of intepolation (p 94), using post #8 as an example

Solve x, such that sinc(x) = k = 1/6.5 ~ 0.153846154

Using my asympotic formula, x ~ Pi/(1 + k + 1.244 k^3) ~ 2.712066499

sinc(2.712) = 0.1535768820, error ~ -0.00027, thus x too high (sinc(x) is decreasing for x = [0, 4.4934])
sinc(2.711) = 0.1539688052, error ~ +0.00013, thus x too low
sinc(2.7115) = 0.1537728267, error ~ -0.00007, used for 3 point quadratic fit.

Acton variation is to sort the values, interpolated estimate up on top.
This ensures top numbers likely the best estimate, with other entries very "tight"
Mistakes during manual calculations are easier to spot.

All entries below are linear interpolation (against top entry), for k = 1/6.5

PHP Code:
sinc(x)       x         0.1537728267  2.71 150.1539688052  2.71 10  1312 9200.1535768820  2.71 20  1312 888  910 --> x = 2.711312910

Quadratic fit might not be needed, if interpolated estimate used as 3rd point
PHP Code:
sinc(x)       x         0.1539688052  2.7110.1535768820  2.712        2.7113129470.1538461395  2.711312947  2.711312910 --> x = 2.711312910

For Casio FX-115MS, Interpolation (A,B), (C,D) Formula:
With CALC: 0A + B + (C-A)-1(D-B)(X-A)
Or, SOLVE: (X-A)(B-0C-D) + (Y-B)(C-A)

update: instead of manually getting the points to fit, we can use iterations:

sinc(pi/2) = 1/(pi/2) = 2/pi ≈ 0.6366197724
If k ≥ 2/pi: x = sin(x)/k
If k < 2/pi: x = pi - asin(kx)

For above example, k < 2/pi, we get 2.712 → 2.7111966 → 2.711332599 ...
Using Aitken Δ² method, above 3 numbers extrapolated to 2.711312910
10-01-2019, 06:03 PM (This post was last modified: 03-19-2021 08:37 PM by Albert Chan.)
Post: #13
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
Inspired by Happy Coincidence for the approximate solution of x tan(x) = k.

I tried the same trick for sin(x)/x = k, using Mathematica.

> <<CalculusPade
> sinc[x_] := Sin[x] / x
⇒ (166320 - 22260 x^2 + 551 x^4) / (166320 + 5460 x^2 + 75 x^4)

> Solve[166320 - 22260 x^2 + 551 x^4 == 0, x] // N
⇒ {{x → -3.14572}, {x → 3.14572}, {x → -5.52302}, {x → 5.52302}}

$$\Large { x^2-\pi^2 \over sinc(x) } = { x^2-\pi^2 \over k}$$

We don't want expression too complicated, and limit it to solving Quadratics:

> lhs = Pade[(x^2-Pi^2) / sinc[x], {x,0,4,2}];
> soln = Solve[ Numerator[lhs] * k - Denominator[lhs] * (x^2-Pi^2) == 0, x]

Although it look like a quartic, it is really a quadratic, with variables x^2.

Instead of solving asinc(k) for x, I choose solving asinc(1-y) for x.
This is a rule I learned from Acton's book, "Numerical Methods that works", p71
Quote:Store the small quantites; Compute the larger ones; Never subtract nearly equal quantities.

$$\large asinc(1-y) ≈ \Large \sqrt{\frac{12y}{1-0.20409y \;+ \sqrt{1 - 0.792y - 0.0318y^2}}}$$

The estimate formula is good for full y range.
y = 1.00 → worst over-estimated: abs error = 0.00023, rel error = 0.0075%
y = 0.85 → worst under-estimated: abs error = 0.00021, rel error = 0.0075%

Previous post example is close to worst case, y = 1-1/6.5 ≈ 0.846154
Using above formula, we get asinc(1/6.5) ≈ 2.71111, under-estimated 0.00020
06-18-2020, 11:54 PM
Post: #14
 Albert Chan Senior Member Posts: 2,066 Joined: Jul 2018
RE: arcsinc( 1-y ), for small y
Re-phrase previous post asinc(1-y) in terms of k, with shorter constants:

$$\large asinc(k) ≈ \Large \sqrt{{12(1-k) \over 1-0.2041(1-k)\;+ \sqrt{(27.11-k)(0.0065+0.0318k)}}}$$

Free-42 code, Newton's method using above guess:

Code:
00 { 85-Byte Prgm } 01▸LBL "ASINC" 02 27.11 03 X<>Y 04 - 05 LASTX 06 STO ST Z     ; k   27.11-k   k 07 0.0318 08 × 09 65ᴇ-4 10 + 11 × 12 SQRT 13 1 14 STO- ST Z 15 +            ; 1+sqrt1   k-1 k 16 0.2041 17 RCL× ST Z 18 + 19 ÷ 20 -12 21 × 22 SQRT 23▸LBL 01       ; newton, f = sinx - kx 24 RCL× ST Z    ; kx    .   k   . 25 LASTX 26 STO ST Z     ; x     kx  x   k 27 SIN 28 - 29 RCL ST Y     ; x    -f   x   k 30 COS 31 RCL- ST T    ; f'   -f   x   k 32 ÷            ; e    x    k   k 33 STO+ ST Y    ; e    x'   k   k 34 X↑2 35 X<>Y 36 STO+ ST Y    ; x' x'+e^2 k   k 37 X≠Y? 38 GTO 01 39 END

Redo previous posts examples:

asinc(5280/5281) → 0.03370775880987942934821849830508089
asinc(0.99) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿→ 0.2453178088540252967075843126095618
asinc(5280/5380) → 0.3348901065998630063407861032211992
asinc(1/6.5) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿→ 2.711312910356983336479873410532187
 « Next Oldest | Next Newest »

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