# How is optimal coordinates change chosen for Chebyshev expansion?

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:

$$c_m=\frac2\pi\int_{-1}^1\frac{T_m(x)J_0\left(\sqrt{8(x+1)}\right)}{\sqrt{1-x^2}}dx,$$

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?

• 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)$$ – uranix Jun 6 '15 at 11:42
• @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. – Ruslan Jun 6 '15 at 11:45
• 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$. – 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.