I'm looking into SLATEC implementation of Bessel function $J_0$ computation (readable in C in GSL), namely at its part for arguments in interval $(0,4)$. There a Chebyshev expansion is used, but the input variable is transformed in non-obvious way. As I've understood, their coefficients are computed using this formula:


and this change of coordinates from $x$ to $\sqrt{8(x+1)}$ is meant to transform from working range of Chebyshev polynomials $(-1,1)$ to the interesting range of Bessel function arguments $(0,4)$. But this change seems somewhat complicated. Naively, I'd rather start from something like $2(x+1)$, which also maps $(-1,1)$ to $(0,4)$. But when I tried it, I found that it gives much worse precision with $13$ terms of the expansion (max absolute error $6\times10^{-11}$ vs $3\times10^{-16}$ in SLATEC).

So I wonder now, how was that transformation, which is actually used in the implementation, found? Can it be proved to be optimal for the precision wanted?

  • $\begingroup$ Do I understand correctly that $c_m$ are coefficients of Chebyshev polynomial fit of $J_0(\sqrt{8(x+1)})$? I.e. $$J_0(\sqrt{8(x+1)}) \approx \sum_{m=0}^{12} c_m T_m(x)$$ $\endgroup$ – uranix Jun 6 '15 at 11:42
  • 1
    $\begingroup$ @uranix Right, but for $m=0$ the sum should have $\frac{c_0}2$, not $c_0$ according to the orthogonality relation of Chebyshev polynomials. $\endgroup$ – Ruslan Jun 6 '15 at 11:45
  • $\begingroup$ Just an idea. Minimize absolute value of $c_{13} = \frac{2}{\pi}\int_{-1}^1 \frac{T_{13}(x) J_0(s(x))}{\sqrt{1-x^2}}dx$ with respect to $s(x): s(-1) = 0, s(1) = 4$. $\endgroup$ – uranix Jun 6 '15 at 11:48

Let $s(x)$ be the substitution, i.e. we are looking for an approximation $$ J_0(s(x)) \approx P(x) = \frac{c_0}{2}T_0(x) + \sum_{m=1}^{M} c_m T_m(x), \quad s(-1) = 0,\, s(1) = 4. $$ with $$ c_m = \frac{2}{\pi}\int_{-1}^1 \frac{T_m(x) J_0(s(x))}{\sqrt{1-x^2}} dx. $$

Let's take such $s(x)$ so $J_0(s(x))$ is a polynomial $Q(x)$ of degree $M$ or less. Now $J_0(s(x))$ is approximated by $P(x)$ exactly. That $s(x)$ would be $$ s(x) = J_0^{-1}(Q(x)). $$ Now the problem is not about optimal $s(x)$ (which there are infinitely many) but about simple $s(x)$ of that form. It should not at least involve any special functions.

By taking $$ J_0(x) = 1 - \frac{x^2}{4} + O(x^4) $$ and neglecting $O(x^4)$ term we can approximate $J_0^{-1}(v)$ as $$ J_0^{-1}(v) \approx 2\sqrt{1-v}. $$ So $$ s(x) = 2\sqrt{1-Q(x)}\\ s(-1) = 2 \sqrt{1 - Q(-1)} = 0\\ s(1) = 2 \sqrt{1 - Q(1)} = 4\\ Q(-1) = 1, \quad Q(1) = -3. $$ A linear $Q(x)$ would be $$ Q(x) = -1 - 2x. $$ so $$ s(x) = 2 \sqrt{1 + 1 + 2x} = 2 \sqrt{2 + 2x}. $$ That's exactly the substitution from GSL.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.