[math-fun] Oscillation and the perils of numerical pessimism

7 May 2012

      The Minsky recurrence handily generates long sequences of a*sin(n*ω+φ).
One can initialize the general case

minsksin[ω_, φ_: 0, ε_: 2 Sin[ω/2], a_: 1] :=
 {a*Sin[φ], -(a/2) ε Cos[φ + ω/2] Csc[ω/2], (4 Sin[ω/2]^2)/ε, ε}

where ε is arbitrary.  E.g.,

In[138]:= minsksin[ω, φ, ε, a]

Out[138]= {a Sin[φ], -(1/2) a ε Cos[φ + ω/2] Csc[ω/2], (4 Sin[ω/2]^2)/ε, ε}

The general Minsky step is
stepminsk[{x_, y_, δ_, ε_}] := {#, y + ε*#, δ, ε} &[x - δ*y]

E.g., running three steps,

In[139]:= Simplify[Nest[stepminsk, %, 3]]

Out[139]= {a Sin[φ + 3 ω], -(1/2) a ε Cos[φ + (7 ω)/2] Csc[ω/2],
 4 Sin[ω/2]^2/ε, ε}

(The 3rd and 4th positions are just along for the ride.)

Nice choices of ε are 2 Sin[ω/2] and 1, where the latter saves a multiply
in the stepper:
stepm1nsk[{x_, y_, δ_}] := {#, y + #, δ}&[x - δ*y]

E.g., for amplitude 1, phase 0,
In[140]:= Drop[minsksin[w, 0, 1], -1]

Out[140]= {0, (-(1/2))*Cot[w/2], 4*Sin[w/2]^2}

Run four steps:
In[144]:= FullSimplify[Nest[stepm1nsk, Out[140], 4]]

Out[144]= {Sin[4*w], (-(1/2))*Cos[(9*w)/2]*Csc[w/2], 2 - 2*Cos[w]}

But instead of Minsky, we could use good old Chebychyov*:

In[37]:= chebsin[ω_, φ_: 0] := {Sin[φ], Sin[φ - ω], 2*Cos[ω]}

whose stepper is

stepcheb[{x_, y_, δ_}] := {δ*x - y, x, δ}

E.g.,
In[146]:= FullSimplify[Nest[stepcheb, chebsin[a], 3]]

Out[146]= {Sin[3 a], Sin[2 a], 2 Cos[a]}

So let's run cheb vs minsky a few million and see who's better.
Using ω := π φ (which will figure in a later email):

In[148]:= chebsin[N[π*GoldenRatio]]

Out[148]= {0, 0.932032, 0.72475}

In[149]:= Timing[Nest[stepcheb, %, 9999999]]

Out[149]= {22.9417, {0.748998, 0.888957, 0.72475}}

In[150]:= minsksin[N[π*GoldenRatio], 0, 1]

Out[150]= {0, 0.730862, 1.27525, 1}

In[151]:= Timing[Nest[stepm1nsk, Drop[%, -1], 9999999]]

Out[151]= {40.471, {0.748998, 0.858748, 1.27525}}

Well, cheb is faster, but which is better?

In[152]:= Sin[9999999*π*GoldenRatio] - %149[[2, 1]]

Out[152]= 3.01301*10^-9

In[153]:= Sin[9999999*π*GoldenRatio] - %151[[2, 1]]

Out[153]= 3.01292*10^-9

They're the same!  But wait, this is a trap.  You can't accurately compute
sin(huge) with single precision!  sin(huge) is like asking which way
your valve stem will point after you drive across the country.

In[154]:= Sin[9999999*π*GoldenRatio] - SetPrecision[%149[[2, 1]], 22]

Out[154]= 1.9186385802598*10^-9

Significantly different, but

In[155]:= Sin[9999999*π*GoldenRatio] - SetPrecision[%151[[2, 1]], 22]

Out[155]= 1.9185395483660*10^-9

They're *still* the same!  Something else is fishy.  MachinePrecision is
~16,
so these values are ~ 10^6 LSBs, not 10^(6/2).  Let's drive ten times
further.

In[156]:= Timing[Nest[stepcheb, %149[[2]], 90000000]]

Out[156]= {214.079, {-0.722388, -0.906264, 0.72475}}

and now we're off by

In[159]:= Sin[99999999*π*GoldenRatio] - SetPrecision[%156[[2, 1]], 22]

Out[159]= -2.00240499632334*10^-8

Ten times as far!  Whereas with Minsky

In[157]:= Timing[Nest[stepm1nsk, %151[[2]], 90000000]]

Out[157]= {382.448, {-0.722388, -0.866577, 1.27525}}

In[160]:= Sin[99999999*π*GoldenRatio] - SetPrecision[%157[[2, 1]], 22]

Out[160]= -2.00242976539902*10^-8

The same!

The explanation:  These recurrences don't compute Sin[n*π*GoldenRatio],
but rather Sin[n*N[π*GoldenRatio]], which is off by

In[161]:= SetPrecision[N[π*GoldenRatio], 22] - π*GoldenRatio

Out[161]= -3.34860*10^-16

We are doomed to linearly creeping error.  But hey, we can tack on a dozen
extra digits
and promise not to drive past a trillion.  But first, let's drive just 999
to see how
long it will take.  (This is where it hits the fan.)

In[166]:= Nest[stepcheb, chebsin[N[π*GoldenRatio, 33]], 999]

Out[166]= {0.*10^121, 0.*10^121, 0.72474978016096023991729327494997}

WAH?

In[167]:= Nest[stepm1nsk, Drop[minsksin[N[π*GoldenRatio, 33], 0, 1], -1],
999]

Out[167]= {0.*10^434, 0.*10^434, 1.27525021983903976008270672505003}

WHA??  This is too gross:

In[168]:= Nest[stepm1nsk, Drop[minsksin[N[π*GoldenRatio, 288], 0, 1], -1],
999]

Out[168]= {0.*10^179, 0.*10^179,
1.27525021983903976008270672505002620127826891198880290709968642251975\
2987505102052895607541471314178683631575021864167072493028631467693604\
7459492180699315909580232204640529914388966962377507409450195194836649\
3221557106208837998002035152821718325478473977678977255669689503145798\
805151318}

"Precision Tracking" (precision trashing, interval arithmetic) *totally*
destroys the computation.  Whereas the actual loss during the fixed
machine precision iteration was negligible.
--rwg

*(The standard counterexample to the Continuum Hypothesis:  The cardinal
between
 ℵ0 and C is the number of ways to spell Chebychyov.)

Bill Gosper

Bill Gosper

Bill Gosper

tags

participants (1)