# A series for the inverse sine cardinal function

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

1. ### David W. CantrellGuest

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

2. ### G. A. EdgarGuest

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

3. ### Rob JohnsonGuest

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
4. ### Peter O. BrackettGuest

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