30
$\begingroup$

Bug introduced in 3.0 and persisting through 12.0 (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?

$\endgroup$
  • 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 '15 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 '15 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 '15 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 '15 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 at 15:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.