Thanks to @steppi for the great work integrating the angular Mathieus into SciPy. Now I’d like to raise a second, smaller break for the Mathieu family while the 2.0 window is open: making the angular argument radians rather than degrees, and fixing the derivative convention at the same time. My original code used radians. My case:
Degrees are an artifact of the old Zhang & Jin implementation (or maybe the old SciPy port). The new kernels are radian-native; the degree convention now exists only as a conversion applied at the SciPy boundary. Nothing downstream of the API wants it.
mathieu_cem and mathieu_sem are the only angular-argument functions in scipy.special that take degrees. Every other trigonometric-flavored function in the module takes radians, as do the corresponding functions in Mathematica, GSL, and Boost.
There is a correctness problem entangled with this. Issue #14580 (open since 2021) reports that the returned derivative is d/dx with respect to the radian angle, while the function value is evaluated at a degree argument. Users must apply a factor of pi/180 by hand, and nothing in the docstring says so. Anyone feeding these into an optimizer or a sensitivity analysis silently picks up a factor of 57.3.
Regarding #4, I’d argue for fixing both at once. If we change only the derivative, that’s a silent numerical break for everyone who has already worked around #14580. If we change only the units, we leave the mismatch in place. Changing both together in the 2.0.0 release is one big break that any Mathieu user cannot fail to notice and fix. And a change in the revision level is the right time to make breaking changes which improve the codebase.
Yes, I agree. Thanks for bringing this back up.The peculiar mixed behavior for the derivatives is especially pernicious.
Some incidentals that I just want to mention for the sake of posterity. For 1) I added both a radian native and degree native path to the new kernels, making them switchable with a template argument. At the SciPy boundary, the degrees path is selected.
For 3), I have a PR that I will submit shortly which will undo this breaking behavior. The wrappers will behave like ufuncs again.
ii) Add a keyword argument for selecting the angle units. Start with default “legacy” (we can’t really call it degrees, because the behavior of the derivative is unusual), and warn if it is selected. Change the default to “radians”, and then deprecate the kwarg.
ii) is a headache, so I am partial to i). I wouldn’t want to make this change directly even in a new major release.
If we come up with a naming scheme, it should be consistent across all of the angular and radial functions.
Regarding which of the two alternatives to take, that’s something for the SciPy team to decide. I think that has more to do with release policy questions than math. But if you go with renaming (option i) I suggest naming the new Mathieu fcns as follows:
# Eigenvalues:
mathieu_a(m, q) # No need to change these since degrees are not involved.
mathieu_b(m, q)
# Angular fcns:
mathieu_ce(m, q, z) # ce is the actual name of this fcn (“cosine elliptic”) per the DLMF and other references.
mathieu_se(m, q, z) # “sine elliptic”
# I never understood the name mathieu_cem. I always figured the “m” in cem meant “Mathieu” but it’s
# redundant since the fcns are already called mathieu_cem with mathieu prefix.
# Reference to DLMF names: https://dlmf.nist.gov/28.2#vi
# Radial fcns:
mathieu_mc1(m, q, z) # Even radial fcn of first type
mathieu_ms1(m, q, z) # Odd radial fcn of first type
mathieu_mc2(m, q, z) # Even radial fcn of second type
mathieu_ms2(m, q, z) # Odd radial fcn of second type
# This follows the DLMF convention: https://dlmf.nist.gov/28.20#iv
Note that some sources (e.g. Wikipedia) use Ce and Se for the radial fcns. However, the DLMF uses Mc and Ms and I take it as authoritative. I communicated with the editor of the Mathieu page
on the DLMF and he said Mc and Ms are the actual fcns. Ce and Se have a different normalization. I tried to follow the DLMF in all the work I did on these things … the nomenclature used in the literature is all over the place but I think the DLMF should be used as the standard.
I like your suggested names. I believe mathieu_cem comes from the notation ce\_{m} for an even Mathieu function of order m and similarly for the others.