Repository: scipy/scipy
Pull request: DOC: `special.mathieu_c/sem`: add information to docstring
#### Reference issue
Closes gh-4096

#### What does this implement/fix?
Adds requested information and examples to the documentation of `mathieu_cem` and `mathieu_sem`.

*First push with `[docs only]` in the commit message did not build the docs, so force-pushing without that...*

Problem statement:
DOC: special: Clarify Mathieu function documentation
Hi, I came across the Mathieu functions and could not get the values from scipy.special.mathieu_cem in accordance with the solutions of the corresponding differential equation.

for q == 0 one should get simply cos(sqrt(a) \* t), but what scipy.special.mathieu_cem returns is more like cos(a \* t). So I assumed that the order parameter a is just defined differently. but then still the solutions for q != 0 do not math. wolfram alpha gives results close to the ones I got from solving the differential equation.

So I conclude the Mathieu functions in scipy are either wrong or do not correspond to the standard definition (see wikipedia, wolfram)

Cheers, Richard

Hints from issue discussion:
ok I get it now! it helps to look into the book Computation of special functions by Shanjie Zhang, Jianming Jin.

I was mistaken that 'a' as the characteristic value, corresponds to the order m used by mathieu_cem. It is rather that way, that for a given parameter 'q' we can calculate the characteristic value 'a' such that the solutions become periodic with period 2pi. And with that condition fulfilled, there exists a hole bunch of solutions, named by the integer parameter m.

So I guess some additional information in the doc of scipy.special.mathieu_cem and as well scipy.special.mathieu_sem would be helpful to know how to use these functions.

Richard

The current version of the documentation strings is here in the sources:
https://github.com/scipy/scipy/blob/master/scipy/special/add_newdocs.py#L1921

If you want to suggest an addition to include in the documentation,
please write it down here. I don't immediately recall the convention
used for these functions here, so if you know how to clarify it that
would be useful.

I think an addition would help a lot. so I would suggest to add here
https://github.com/scipy/scipy/blob/master/scipy/special/add_newdocs.py#L1941

```
    Even Mathieu function and its derivative

    Returns the even Mathieu function, ``ce_m(x,q)``, of order m and
    parameter q evaluated at x (given in degrees).  Also returns the
    derivative with respect to x of ce_m(x,q)

    Notes
    -----
    The Mathieu function obeys the following differential equation:

    .. math:: \mathrm{ce}_\mathrm{m}''(x) + [\alpha + 2 q \cos(2x)] \mathrm{ce}_\mathrm{m}(x) = 0

    For a given real parameter 'q' the characteristic value 'alpha' is chosen
    such that ce_m(x) is periodic with period 2 pi (see mathieu_a). The possible
    solutions for the characteristic value 'alpha' are indexed with 'm'. Further the even
    function fulfills the initial condition ce_m(x=0) != 0 and ce_m'(x=0) == 0 where as
    for the odd function ce_m(x=0) == 0 and ce_m'(x=0) != 0 holds true.

    Further reading see [1] page 475.

    References
    ----------
    .. [1] Shanjie Zhang and Jianming Jin, Computation of special functions, John Wiley, 1996
```

There seems to be still some contradictions, because the reference claims for the even functions to have a period of pi, where a simple plot shows that the scipy mathieu functions have 2pi period, both, even and odd. Further the value for ce_m(x=0) is actually arbitrary as we have a homogeneous ODE. So where does the actual value for that implementation come from?

The addition could also be added here
https://github.com/scipy/scipy/blob/master/scipy/special/add_newdocs.py#L2040

ok enough for now ...

best, Richard

The Mathieu function conventions are the same as in Abramowitz&Stegun, if I remember correctly from checking this a couple of years ago. There is also a test checking that the values produced are in accordance with the implementation in GSL.

Will probably require some research, but once the math is understood, doesn't seem like it would be hard to address in the documentation.
I can see how correspondence with the Mathworld reference is not as obvious, but I found the function to correspond closely with information on Wikipedia.

> for q == 0 one should get simply cos(sqrt(a) * t),

Wikipedia states that they should be proportional to $\cos(nx)$ and $\sin(nx)$.

<img width="1388" alt="image" src="https://github.com/user-attachments/assets/ee1846cd-0e83-4feb-998b-ea33b917d55c">

Indeed, they are, if we recognize `m` to correspond with $n$.

<details>

```python3
import numpy as np
from scipy import special
import matplotlib.pyplot as plt

m = 4
q = 0
x = np.linspace(0, 360, 300)

y, dy = special.mathieu_cem(m, q, x)
plt.plot(x, y)
plt.plot(x, np.cos(m*x/180*np.pi), '--')

y, dy = special.mathieu_sem(m, q, x)
plt.plot(x, y)
plt.plot(x, np.sin(m*x/180*np.pi), '--')

plt.legend(('mathieu_cem', 'cos(mx)', 'mathieu_sem', 'sin(mx)'))
```

</details>

<img width="925" alt="image" src="https://github.com/user-attachments/assets/3d657876-eb90-4cf9-92a4-94429f34ad21">

---

> I was mistaken that 'a' as the characteristic value, corresponds to the order m used by mathieu_cem. It is rather that way, that for a given parameter 'q' we can calculate the characteristic value 'a' such that the solutions become periodic with period 2pi. And with that condition fulfilled, there exists a hole bunch of solutions, named by the integer parameter m.

This is addressed by the definition of $ce_n$ and $se_n$:

<img width="1386" alt="image" src="https://github.com/user-attachments/assets/ae94f371-d21e-4e83-8b9a-e21921690e3f">

---

> There seems to be still some contradictions, because the reference claims for the even functions to have a period of pi, where a simple plot shows that the scipy mathieu functions have 2pi period, both, even and odd

Actually, a plot confirms the property claimed by the reference.

<details>

```python3
import numpy as np
from scipy import special
import matplotlib.pyplot as plt

q = 100
x = np.linspace(0, 360, 300)
y, dy = special.mathieu_cem(np.asarray([1, 3, 5])[:, np.newaxis], q, x)
plt.plot(x, y.T)
plt.legend(('m = 1', 'm = 3', 'm = 5'))
```

</details>

Period $2\pi$ for odd $m$.
<img width="990" alt="image" src="https://github.com/user-attachments/assets/a8358a3d-2cdc-4334-b825-360825fd2572">

Period $\pi$ for even $m$.
<img width="968" alt="image" src="https://github.com/user-attachments/assets/75118f22-c2d7-49a2-b8ba-9edf43ebd292">

---

> Further the value for ce_m(x=0) is actually arbitrary as we have a homogeneous ODE. So where does the actual value for that implementation come from?

<img width="1369" alt="image" src="https://github.com/user-attachments/assets/4fe977c7-cecf-481f-8559-7a523853327b">

---

@j-bowhay, it looks like you began to review gh-17333. I closed it for a few reasons (inaccurate, incomplete, and inactive), but would you be interested in either writing or reviewing a PR that added some documentation to the `mathieu_cem` and `mathieu_sem` functions to close this issue? To keep the scope manageable, I think it should:

* identify the differential equation the function solves
* identify the relationship between the arguments of the function and the variables in the differential equation
* reference Wikipedia for other properties and the normalization scheme
* show one example, e.g. of the $\pi$ or $2\pi$ periodicity depending on $m$.
I am currently away travelling until next year so unfortunately I have quite limited computer access

Commit message: DOC: special.mathieu_c/sem: add information to docstring

Patch to review:
diff --git a/scipy/special/_special_ufuncs_docs.cpp b/scipy/special/_special_ufuncs_docs.cpp
index fedf9fa87a81..91ccc742650d 100644
--- a/scipy/special/_special_ufuncs_docs.cpp
+++ b/scipy/special/_special_ufuncs_docs.cpp
@@ -4762,6 +4762,42 @@ const char *mathieu_cem_doc = R"(
     --------
     mathieu_a, mathieu_b, mathieu_sem
 
+    Notes
+    -----
+    Even Mathieu functions are the solutions to Mathieu's differential equation
+
+    .. math::
+
+        \frac{d^2y}{dx^2} + (a_m - 2q \cos(2x))y = 0
+
+    for which the characteristic number :math:`a_m` results in an even, periodic
+    solution :math:`y(x)`` with period 180 degrees (for even :math:`m`) or
+    360 degrees (for odd :math:`m`).
+
+    References
+    ----------
+    'Mathieu function'. *Wikipedia*.
+    https://en.wikipedia.org/wiki/Mathieu_function
+
+    Examples
+    --------
+    Plot Mathieu functions of orders ``2`` and ``4``.
+
+    >>> import numpy as np
+    >>> from scipy import special
+    >>> import matplotlib.pyplot as plt
+    >>> m = np.asarray([2, 4])
+    >>> q = 50
+    >>> x = np.linspace(-180, 180, 300)[:, np.newaxis]
+    >>> y, _ = special.mathieu_cem(m, q, x)
+    >>> plt.plot(x, y)
+    >>> plt.xlabel('x (degrees)')
+    >>> plt.ylabel('y')
+    >>> plt.legend(('m = 2', 'm = 4'))
+
+    Note that the functions are even. Also, because the orders ``2`` and
+    ``4`` are even, the period of each function is 180 degrees.
+
     )";
 
 const char *mathieu_modcem1_doc = R"(
@@ -4927,6 +4963,42 @@ const char *mathieu_sem_doc = R"(
     --------
     mathieu_a, mathieu_b, mathieu_cem
 
+    Notes
+    -----
+    Odd Mathieu functions are the solutions to Mathieu's differential equation
+
+    .. math::
+
+        \frac{d^2y}{dx^2} + (a_m - 2q \cos(2x))y = 0
+
+    for which the characteristic number :math:`a_m` results in an odd, periodic
+    solution :math:`y(x)`` with period 180 degrees (for even :math:`m`) or
+    360 degrees (for odd :math:`m`).
+
+    References
+    ----------
+    'Mathieu function'. *Wikipedia*.
+    https://en.wikipedia.org/wiki/Mathieu_function
+
+    Examples
+    --------
+    Plot Mathieu functions of orders ``2`` and ``4``.
+
+    >>> import numpy as np
+    >>> from scipy import special
+    >>> import matplotlib.pyplot as plt
+    >>> m = np.asarray([2, 4])
+    >>> q = 50
+    >>> x = np.linspace(-180, 180, 300)[:, np.newaxis]
+    >>> y, _ = special.mathieu_sem(m, q, x)
+    >>> plt.plot(x, y)
+    >>> plt.xlabel('x (degrees)')
+    >>> plt.ylabel('y')
+    >>> plt.legend(('m = 2', 'm = 4'))
+
+    Note that the functions are even. Also, because the orders ``2`` and
+    ``4`` are even, the period of each function is 180 degrees.
+
     )";
 
 const char *modfresnelm_doc = R"(


Code review comments:
diff_hunk: @@ -4927,6 +4963,42 @@ const char *mathieu_sem_doc = R"(
     --------
     mathieu_a, mathieu_b, mathieu_cem
 
+    Notes
+    -----
+    Odd Mathieu functions are the solutions to Mathieu's differential equation
+
+    .. math::
+
+        \frac{d^2y}{dx^2} + (a_m - 2q \cos(2x))y = 0
path: scipy/special/_special_ufuncs_docs.cpp
text: @author:
```suggestion
        \frac{d^2y}{dx^2} + (b_m - 2q \cos(2x))y = 0
```
Maybe do this to be suggestive of `mathieu_b`? 

@user1:
Feels like that might be confusing; we're not solving a different equation.

@author:
Well, it's the same form of the equation but with a different infinite set of parameters that we call $b_m$ instead of $a_m$ (apparently), according to the name of the functions (`mathieu_b` vs `mathieu_a`). Different symbols are sometimes used for different numbers. 

But I would only do this if I also mentioned the function `mathieu_b` below. If I were to include it, it seemed more clear to write that the characteristic parameters $b_m$ (rather than $a_m$) can be calculated with `mathieu_b`.

Perhaps it would have been most clear of the functions had been named differently, e.g. `mathieu_cem_a`/`mathieu_sem_a` or one function with a `ce_m`/`se_m` or "even"/"odd" option, but we are stuck with this letter `b` in the name of one.

diff_hunk: @@ -4762,6 +4762,42 @@ const char *mathieu_cem_doc = R"(
     --------
     mathieu_a, mathieu_b, mathieu_sem
 
+    Notes
+    -----
+    Even Mathieu functions are the solutions to Mathieu's differential equation
+
+    .. math::
+
+        \frac{d^2y}{dx^2} + (a_m - 2q \cos(2x))y = 0
+
+    for which the characteristic number :math:`a_m` results in an even, periodic
+    solution :math:`y(x)`` with period 180 degrees (for even :math:`m`) or
path: scipy/special/_special_ufuncs_docs.cpp
text: @author:
```suggestion
    solution :math:`y(x)` with period 180 degrees (for even :math:`m`) or
```

diff_hunk: @@ -4927,6 +4963,42 @@ const char *mathieu_sem_doc = R"(
     --------
     mathieu_a, mathieu_b, mathieu_cem
 
+    Notes
+    -----
+    Odd Mathieu functions are the solutions to Mathieu's differential equation
+
+    .. math::
+
+        \frac{d^2y}{dx^2} + (a_m - 2q \cos(2x))y = 0
+
+    for which the characteristic number :math:`a_m` results in an odd, periodic
+    solution :math:`y(x)`` with period 180 degrees (for even :math:`m`) or
path: scipy/special/_special_ufuncs_docs.cpp
text: @author:
```suggestion
    solution :math:`y(x)` with period 180 degrees (for even :math:`m`) or
```

diff_hunk: @@ -4762,6 +4762,42 @@ const char *mathieu_cem_doc = R"(
     --------
     mathieu_a, mathieu_b, mathieu_sem
 
+    Notes
+    -----
+    Even Mathieu functions are the solutions to Mathieu's differential equation
+
+    .. math::
+
+        \frac{d^2y}{dx^2} + (a_m - 2q \cos(2x))y = 0
+
+    for which the characteristic number :math:`a_m` results in an even, periodic
+    solution :math:`y(x)`` with period 180 degrees (for even :math:`m`) or
+    360 degrees (for odd :math:`m`).
+
+    References
+    ----------
+    'Mathieu function'. *Wikipedia*.
path: scipy/special/_special_ufuncs_docs.cpp
text: @user1:
Is this valid reference format? I expect:

```
.. [1] O. McNoleg, ...
```

But perhaps SciPy follows its own standard.

@author:
~~Is McNoleg the author of this article?~~(I figured it was a fake name.) There is plenty of precedent in SciPy for not listing an author for Wikipedia. It would typically be like "The Wikipedia contributors" if it were there, which doesn't really add identifying information, so it's often (always? I just searched... looks like always) omitted.

@user1:
Oh, I meant the ".. [1]" bit; I thought that was necessary for numpydoc to parse the reference.

@author:
Oh! I messed up. I'll fix those.