34

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?

xzczd
  • 65,995
  • 9
  • 163
  • 468
Ruslan
  • 7,152
  • 1
  • 23
  • 52
  • 2
    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. – Histograms May 14 '15 at 15:24
  • 1
    @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. – Ruslan May 14 '15 at 17:50
  • 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. – Histograms May 14 '15 at 18:54
  • 1
    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. – Oleksandr R. May 15 '15 at 14:13
  • @OleksandrR. your suggested command fails for z>=6 because of too small $MaxExtraPrecision (messages are visible when doing ListPlot of Table instead of just Plot), but yeah, it seems a bit more robust than working with machine precision. I'll update the question. – Ruslan May 16 '15 at 07:05
  • Just tested in v3, it's broken from the beginning: https://i.stack.imgur.com/dcrVw.png – xzczd Jun 15 '19 at 15:11
  • 4
    @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. – Ruslan Jun 15 '19 at 15:33

1 Answers1

6

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

Henri Leveque
  • 61
  • 1
  • 2