A series for the inverse sine cardinal function

Discussion in 'Numerical Analysis' started by David W. Cantrell, Feb 14, 2004.

  1. The sine cardinal function,

    ( 1 if x = 0
    sinc(x) = (
    ( sin(x)/x otherwise,

    arises in several applications, and so its inverse should also be of
    interest. This note gives a series expansion for the inverse of sinc(x),
    0 <= x <= x0, where x0 (approx. 4.4934) denotes the positive value of x
    at which sinc(x) reaches its absolute minimum.

    [Besides posting this to sci.math, I have also posted it to
    sci.math.num-analysis and comp.dsp in the hope that some people might be
    able to give references to previous mathematical treatments of the inverse
    sinc function. In particular, I'd be interested to know if the series given
    here is already known.]

    With f(x) = 2*x + 3*x^3/10 + 321*x^5/2800 + 3197*x^7/56000 +
    445617*x^9/13798400 + 1766784699*x^11/89689600000 +
    317184685563*x^13/25113088000000 +
    14328608561991*x^15/1707689984000000 +
    6670995251837391*x^17/1165411287040000000 +
    910588298588385889*x^19/228420612259840000000 + ...,

    it can be shown that the desired inverse, abbreviated as Asinc here, is
    given by

    Asinc(x) = Sqrt(3/2) * f(Sqrt(1 - x))

    -----------------------------------------------------------------------
    A simple example:
    Find the positive root of the equation sin(x) = x/2.

    Rewriting the equation as sinc(x) = 1/2,
    the solution is x = Asinc(1/2) = Sqrt(3/2) * f(Sqrt(1/2)).

    Using just those terms shown above in the series for f, we get
    x = 1.8954905... (For comparison, solving by a different method, a much
    more accurate approximation is x = 1.8954942668788...)
    -----------------------------------------------------------------------

    The radius of convergence of the above series for f is slightly larger than
    1.1 . As expected, when x is near that value, convergence is very slow.

    The series for f was obtained by reversion of series, etc. To be more
    specific, for anyone interested, in Mathematica the series may be obtained
    using

    Simplify[Normal[InverseSeries[1 - Series[Sin[x]/x, {x, 0, n}]]/Sqrt[3/2]
    /. x -> x^2], x > 0]

    where n (even) should be chosen to be 1 more than the highest power of x
    desired in the result.

    Unfortunately, I do not know a nice formula for the coefficients in the
    series for f.

    Thoughtful comments will be appreciated.

    David Cantrell
     
    David W. Cantrell, Feb 14, 2004
    #1
    1. Advertisements

  2. David W. Cantrell

    G. A. Edgar Guest

    here it is in Maple...
    z := 1/2*6^(1/2)*u
    AA := 1/3*sin(1/2*6^(1/2)*u)*6^(1/2)/u = 1-y^2
    AA1 := -1/3*sin(RootOf(-_Z+_Z*y^2+sin(_Z)))*6^(1/2)/(-1+y^2)
    2*y+3/10*y^3+321/2800*y^5+3197/56000*y^7
    +445617/13798400*y^9
    +1766784699/89689600000*y^11+317184685563/25113088000000*y^13
    +14328608561991/1707689984000000*y^15
    +6670995251837391/1165411287040000000*y^17
    +910588298588385889/228420612259840000000*y^19
    +1889106915501879285127263/669318078043783168000000000*y^21
    +122684251268939994619239/60571771768668160000000000*y^23
    +86578199631805319180104483967/58899990867852918784000000000000*y^25
    +36790913563978761395277930686421/
    34161994703354692894720000000000000*y^27
    +295479400033606079171291233070109663/
    371437974410411888261201920000000000000*y^29
    +5429197579977012936689051219262425180781/
    9174517967937173640051687424000000000000000*y^31
    +111912180845235829957717919886518013347855667/
    252655431286439247725093999083520000000000000000*y^33
    +21623582556767547163123245489615552651038372109/
    64865414807824606864932296091238400000000000000000*y^35
    +1686950689722579328034933293949678511289600201851/
    6691379632807169971329857912569856000000000000000000*y^37
    +362964877894310955248925785551248746981943835416147/
    1895485357802467420969439750506151936000000000000000000*y^39
    +241119375769652087142546687133690376324455849576103305306113/
    1651328495419246089263940926464174667092459520000000000000000000*y^41
    +62733744440681157939079200023293588467527441653821853933262223/
    561451688442543670349739914997819386811436236800000000000000000000*y^43
    +90541073261492600342265407374539407149059862146954660576097727343/
    1055529174271982100257511040195900447205500125184000000000000000000000
    *y^45
    +126749063504538759034999649107808731423757688715140874157349767781/
    1919143953221785636831838254901637176737272954880000000000000000000000
    *y^47
    +249693053510060031661928074632553530869877798468407353253082594529911/
    489709179353451023617580383462989683786482353963008000000000000000000000
    0
    *y^49
    +O(y^50)
     
    G. A. Edgar, Feb 14, 2004
    #2
    1. Advertisements

  3. David W. Cantrell

    Rob Johnson Guest

    sinc(x) >= -0.21723362821 (minimum at x = +/-4.49340945791). Therefore,
    sqrt(1-sinc(x)) <= 1.10328311335, hopefully the radius of convergence.
    Another way to do this is using the Taylor series for x = sinc(y). We
    want to find y in terms of x where

    y^2 y^4 y^6
    x = 1 - --- + --- - ---- + ...
    6 120 5040

    Define f(z) by

    z^2 z^3
    f(z) = z - --- + --- - ...
    20 840

    oo
    --- k-1 k 6
    = > (-1) z -------
    --- (2k+1)!
    k=1
    so that
    1
    x = 1 - - f(y^2)
    6
    and
    -1
    y = sqrt(f (6-6x))

    For the inverse of f, Mathematica gives

    Normal[InverseSeries[Series[Sum[-(-x)^k 6/(2k+1)!,{k,1,10}],{x,0,10}],x]]

    2 3 4 5 6
    x 2 x 13 x 4957 x 58007 x
    x + -- + ---- + ----- + --------- + ----------- +
    20 525 37800 145530000 16216200000

    7 8 9
    1748431 x 4058681 x 5313239803 x
    ------------- + -------------- + ------------------- +
    4469590125000 92100645000000 1046241656460000000

    10
    2601229460539 x
    ----------------------
    4365681093774000000000

    Although this doesn't yield a simple formula for the general term of the
    series, at least the first 40 terms all have positive coefficients. I
    think there is a good chance that all the terms may share this property,
    but I need to look into this a bit more.

    Rob Johnson <>
    take out the trash before replying
     
    Rob Johnson, Feb 15, 2004
    #3
  4. David:

    In general a rational polynomial approximation [mini-max] determined using
    say
    Remez second method will be much more economic than any power series.

    I have used rational polynomial approximations in digital signal processing
    (DSP)
    applications of the sinc function and its' inverse to great practical
    advantage in
    terms of demands upon real time storage and execution cycles.
     
    Peter O. Brackett, Feb 17, 2004
    #4
    1. Advertisements

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.