Software 8.5.1 (The Airy function from its series). The recursion of Example 8.1.2, with initial values \(\operatorname {Ai}(0)\) and \(\operatorname {Ai}'(0)\) taken from SciPy, summed to \(60\) terms and compared with SciPy’s airy:

>>> import numpy as np
>>> from math import factorial
>>> from scipy import special
>>> def airy_series(a0, a1, t, K=60):      # y'' = t y, power series at 0
...     a = [a0, a1, 0.0]
...     for k in range(1, K):
...         a.append(a[k - 1] / ((k + 2) * (k + 1)))
...     return sum(c * t**j for j, c in enumerate(a))
...
>>> c1, c2 = special.airy(0)[0], -special.airy(0)[1]   # Ai(0), -Ai'(0)
>>> for t in (-4.0, 0.5, 2.0):
...     print(t, airy_series(c1, -c2, t), special.airy(t)[0])
...
-4.0 -0.07026553294928879 -0.07026553294928964
0.5 0.23169360648083345 0.23169360648083343
2.0 0.034924130423274184 0.03492413042327436

Reading the output. The series agrees with the library value to about \(14\) significant digits at all three points. At \(t = -4\), where the solution oscillates, the sum is made of terms as large as \(4\) in modulus that cancel to \(-0.07\), and the cancellation grows rapidly with \(|t|\). For large \(|t|\) the series, although convergent, becomes useless numerically because of cancellation, and libraries switch to asymptotic expansions there.

MATLAB/Octave: airy(0, [-4 0.5 2]).

Software 8.5.2 (Bessel functions). Partial sums of the series for \(J_0(5)\), the first zeros of \(J_0\), and the closed form of \(J_{1/2}\):

>>> J0 = lambda t, K: sum((-1)**k * (t / 2)**(2 * k) / factorial(k)**2
...                       for k in range(K))
>>> [round(J0(5.0, K), 6) for K in (5, 10, 15, 20)]
[0.383036, -0.177603, -0.177597, -0.177597]
>>> round(special.jv(0, 5.0), 6)
-0.177597
>>> special.jn_zeros(0, 5)
array([ 2.40482556,  5.52007811,  8.65372791, 11.79153444, 14.93091771])
>>> np.diff(special.jn_zeros(0, 5)), np.pi
(array([3.11525255, 3.1336498 , 3.13780653, 3.13938327]), 3.141592653589793)
>>> t = 3.0; special.jv(0.5, t), np.sqrt(2 / (np.pi * t)) * np.sin(t)
(0.06500818287737586, 0.06500818287737577)

Reading the output. Five terms of the series give \(0.383\), nowhere near \(J_0(5) = -0.177597\); ten give it to five decimals, and fifteen to all six shown. The zeros of \(J_0\) are \(2.405, 5.520, 8.654, \dots \), and their spacings \(3.115, 3.134, 3.138, 3.139\) approach \(\pi \), as the asymptotic form predicts. \(J_{1/2}(3)\) equals \(\sqrt {2/(3\pi )}\sin 3\), confirming Theorem 8.3.1(iv).

MATLAB/Octave: besselj(0, 5); the zeros with fzero(@(t) besselj(0, t), 2.4).

Software 8.5.3 (Legendre polynomials). Rodrigues’ formula in SymPy, checked against the built-in legendre, and the matrix of inner products \(\int _{-1}^1P_mP_n\):

>>> import sympy as sp
>>> x = sp.symbols('x')
>>> rod = lambda n: sp.expand(sp.diff((x**2 - 1)**n, x, n)
...                           / (2**n * factorial(n)))   # Rodrigues
>>> [rod(n) for n in range(4)]
[1, x, 3*x**2/2 - 1/2, 5*x**3/2 - 3*x/2]
>>> rod(3) == sp.expand(sp.legendre(3, x))
True
>>> [[sp.integrate(rod(m) * rod(n), (x, -1, 1)) for n in range(4)]
...  for m in range(4)]
[[2, 0, 0, 0], [0, 2/3, 0, 0], [0, 0, 2/5, 0], [0, 0, 0, 2/7]]

Reading the output. Rodrigues’ formula gives \(1\), \(x\), \(\frac 32x^2 - \frac 12\) and \(\frac 52x^3 - \frac 32x\), the first four Legendre polynomials, and agrees with SymPy’s own. The matrix of integrals is diagonal, with entries \(2, \frac 23, \frac 25, \frac 27\), which are \(\frac {2}{2n + 1}\): Theorem 8.4.1(ii) exactly.

MATLAB/Octave: legendreP(3, x) with the Symbolic Math Toolbox.