34
$\begingroup$

Bug introduced in 3.0 and persisting through 13.2 (reported as CASE:3208982)

I'm trying to plot MathieuC[-3,0.3,I x] for $x\in[0,10]$, and here's what I get even with arbitrary precision arithmetic (here I use ListPlot and Table instead of Plot to make computation faster and use a regular grid):

$MaxExtraPrecision=200;
ListPlot[Table[{z,N[Re@MathieuC[-3,3/10,I z],50]},{z,0,7,1/100}],
  PlotRange->{-0.6,0.6}, Joined->True]

bad behavior of MathieuC

So, starting from about z==3.3, the results are unreliable at all. Comparing machine precision and arbitrary precision computations with N[], I get this:

N[Re@MathieuC[-3, 3/10, 5 I]]
N[Re@MathieuC[-3, 3/10, 5 I], 50]

-1.26013246174486*10^28

-1.2601324617438073657004964476674284869363635740962*10^28

So, it seems not even lack of working precision — looks like the algorithm is broken. To make sure I'm not misunderstanding the supposed behavior of Mathieu functions, I've also checked it by solving Mathieu equation numerically:

With[{a = -3, q = 0.3},
 sol = NDSolveValue[{
      -w''[z] - 2 q Cosh[2 z] w[z] == -a w[z],
      w[0] == Re@MathieuC[a, q, 0],
      w'[0] == 0
      }, w, {z, 0, 7}, MaxSteps -> 10^6]
 ];
$MaxExtraPrecision=200;
Show[
 Plot[sol[z],{z,0,7}, PlotStyle->Darker@Green, PlotPoints->300],
 ListPlot[Table[{z,N[Re@MathieuC[-3,3/10,I z],50]},{z,0,7,1/100}],Joined->True],
PlotRange->{-0.6,0.6}]

MathieuC versus NDSolve result

This shows that the function indeed should look nicer, but apparently MathieuC can't calculate it for even moderately large imaginary arguments. It also appears that Mathematica versions 5 to 12 give exactly the same results (sometimes with differences in several least significant digits).

Is it a bug, or is it documented somewhere? Are there any workarounds, allowing me to evaluate Mathieu functions not as slowly as via NDSolve, and still not reimplementing them like e.g. here?

Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ Yes, I think it is broken. But this may be because MathieuC requires exact arguments. Plot feeds it machine precision values and the documentation says it requires exact values under the Properties&Relations section. $\endgroup$
    – Histograms
    May 14, 2015 at 15:24
  • 1
    $\begingroup$ @Histograms where did you see such a requirement? I checked built-in docs in MMA 9 and 10 as well as online documentation and didn't see it. The only thing in P&R section is that for some machine-precision arguments the precision is insufficient, and one has to use higher precision in N (or just supply numbers with higher precision, not necessarily exact). In any case, I've given an example of exact arguments and N in the OP, with the same wrong results. $\endgroup$
    – Ruslan
    May 14, 2015 at 17:50
  • $\begingroup$ That's what I was referring to and requirement was a strong word. Perhaps for these choices of the other arguments, the precision problems are amplified, I'm just guessing. $\endgroup$
    – Histograms
    May 14, 2015 at 18:54
  • 1
    $\begingroup$ It is a good question. Just to note that it is easier and more convincing to use Plot[Re@MathieuC[-3, 3/10, I z], {z, 0, 10}, PlotRange -> {-0.6, 0.6}, WorkingPrecision -> 50] as the comparison, rather than the one calculated at fixed precision. $\endgroup$
    – Oleksandr R.
    May 15, 2015 at 14:13
  • 4
    $\begingroup$ @xzczd what's interesting, when I tried Maple 2015 to see if it has a better implementation (the functions are of the same type there as in MMA: MathieuC and MathieuS), it appeared to fail in exactly the same way. So I suspect Wolfram and Maplesoft used algorithms from the same paper, which appear broken as published, and no one had noticed this. $\endgroup$
    – Ruslan
    Jun 15, 2019 at 15:33

1 Answer 1

Reset to default
6
$\begingroup$

I, you are using imaginary argument, which is suppose that it's rather a radial Mathieu Function instead of an angular one, Radial Mathieu Function are not implemented in Mathematica and are incompletely simulated by using imaginary argument on an angular Mathieu Function. Angular Mathieu Functions is defined on [0,2*Pi] and have some periodic properties... It's not surprising that functions diverge when argument > 6, since probably Mathematica use matrix or continued fraction methods. To have better estimate of the Radial mathieu functions (so called modified Mathieu Function), one uses Series of cross-product of Bessel functions. But it means that you have to deal with more theoretical matter about these functions and not just use incomplete implementation of Mathieu Function in Mathematica.

Henri Lévêque

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.