Problem with scaling of derivative output of Mathieu functions.

[WarrenWeckesser]

Over in scipy/scipy#14580, it was pointed out that the scaling for the derivative output of scipy.special.mathieu_cem and scipy.special.mathieu_sem was not correct. In those Python functions, the x argument is given in degrees. The problem is that the returned derivative behaves as if the input was in radians--it is missing the scaling factor of π/180.

I have created this issue here because the problem was brought over to xsf. This program (using C++20)

#include <iomanip>
#include <iostream>
#include <numbers>
#include <xsf/mathieu.h>

int main(int argc, char *argv[])
{
    using std::numbers::pi_v;
    double f, df;

    double m = 3.0;
    double q = 1.5;
    double theta = 1.0;                     // radians
    double angle = theta*180/pi_v<double>;  // degrees
    xsf::cem(m, q, angle, f, df);

    std::cout << std::setprecision(15);
    std::cout << "f = " << f << std::endl;
    std::cout << "df          = " << df << std::endl;

    double h = 1e-8;
    double fplus, fminus, tmp;
    xsf::cem(m, q, angle + h, fplus, tmp);
    xsf::cem(m, q, angle - h, fminus, tmp);
    double approx_diff = (fplus - fminus)/(2*h);
    std::cout << "finite diff = " << approx_diff << std::endl;
    std::cout << "df*π/180    = " << df * pi_v<double> / 180.0 << std::endl;

    return 0;
}

prints

f = -0.867785570960645
df          = -1.05398845039289
finite diff = -0.0183955628507704
df*π/180    = -0.0183955687373489

You can check with Wolfram Alpha (MathieuCPrime(MathieuCharacteristicA(3, 3/2), 3/2, 1)) that df is the value that would be correct if the implementation of xsf::cem accepted radians instead of degrees. With degrees, the values of df should be -0.018395568737348884.

The program computes a numerical approximation to the derivative to verify that df is off by a factor of π/180.

The straightforward fix is to multiply the derivative outputs by pi/180 before returning.

(Personally, I wouldn't mind if we changed all the Mathieu functions to accept radians instead of degrees. Over in SciPy, we would have to find a way to make the transition that is not too painful--may just create a new set of functions, and slowly deprecate (or at least doc-deprecate) the ones that accept degrees.)

[brorson]

I'll chime in here and agree whole-hartedly -- you are right. Although I have not appeared her for a couple of months I have some students who have been working on testing the Mathieu implementations -- both the existing buggy Scipy one and the new one in pull request #60. They also found that the Mathieu derivatives are incorrectly scaled. The call assumes the inputs are degrees. The Mathieu return is fine since it is not scaled, but the Mathieu derivative is scaled as if the call was made in radians.

I also cast my vote to change all Mathieu inputs to radians. Those are the natural units. I can't imagine anybody sophisticated enough to use Mathieu functions would want to use inputs in degrees.

Once the semester is over I will have more time to work on the Mathieu project again.

[brorson]

By the way, the impl I offered assumes everything is in rads, so doesn't have this problem.