Mathematica doesn't evaluate MathieuC function for large arguments

Question

Various MathieuC problems identified here have been reported as CASE:3698818.

I need to calculate some numeric integrations of the periodic Mathieu's Functions. So I used NIntegrate and it works well as long as the arguments of the function are not too large. Whenever I try to do the calculation with a large argument, Mathematica leaves the function unevaluated and I get the "The integrand has evaluated to non-numerical values" error.

I noticed that when I round off the arguments, or use Rationalize to get a rational approximation of them, it does evaluate the function, but the problem with that is that the functions are very sensitive to argument changes - as soon as I round off the argument the function stops being periodic and becomes exponentially growing - which is definitely not what I'm looking for. I tried increasing the precision of the approximation, and it wasn't precise enough when Mathematica stopped evaluating the function again.

So what I need is to keep the precision but make Mathematica evaluate the function. (The built-in function Evaluate doesn't work either)

Is there a way to do that?

Here are some code lines for example:

    b = 0.49
    q = q /. FindRoot[-2 q b == MathieuCharacteristicA[6, q], {q, 1000}]
    q2 = q /. FindRoot[-2 q b == MathieuCharacteristicA[8, q], {q, 1000}]
    NIntegrate[MathieuC[-2 q b, q, x]*MathieuC[-2 q2 b, q2, x], {x, 0, 2 Pi}]
    Plot[Evaluate[{Re[#], Im[#]} &@MathieuC[-2 q2 b, q2, x]], {x, 0, 4 Pi}, PlotStyle -> {Blue, {Blue, Dashed}}]

Thanks!

Answer 1

It wasn't real clear to me exactly what the OP did about precision, but code below looks like it might be working. You seem to need fairly high precision for all arguments.

b = 0.49`100;
q1 = q /. 
  FindRoot[-2 q b == MathieuCharacteristicA[6, q], {q, 1000}, 
   WorkingPrecision -> 100]
q2 = q /. 
  FindRoot[-2 q b == MathieuCharacteristicA[8, q], {q, 1000}, 
   WorkingPrecision -> 100]
(*
605.789297872353700695280254027183815298851997479216300967816087036777\
2075840329756326461722267299774

1036.13835973661722165073772712372558794514088740162366560455284235128\
4742499772811507801258473196255
*)

NIntegrate[
 MathieuC[-2 q1 b, q1, x]*MathieuC[-2 q2 b, q2, x], {x, 0, 2 Pi}, 
 PrecisionGoal -> 8, WorkingPrecision -> 100]
(*
1.43071586725038018235424085351766966962166226308387215510919976577346\
8232618999353484020990207944397
*)

data = Table[MathieuC[-2 q1 b, q1, x]*MathieuC[-2 q2 b, q2, x], {x, 0, 2 Pi, 2*Pi/1000}];
plot = ListLinePlot[data, PlotRange -> All, DataRange -> {0, 2 Pi}]

To use Plot instead of ListLinePlot, set WorkingPrecision -> 100. One should also limit MaxRecursion because computation of the Mathieu functions at high precision is slow. The setting MaxRecursion -> 2 reduced the time by about 1/10 of the default setting, from 86 sec. to about 8.7. (PlotPoints -> 100 was used to capture all the oscillations.)

Plot[MathieuC[-2 q1 b, q1, x]*MathieuC[-2 q2 b, q2, x], {x, 0, 2 Pi},
  PlotPoints -> 100, MaxRecursion -> 2,   (* for speed *)
  WorkingPrecision -> 100, PlotRange -> All] // AbsoluteTiming

At least we can show that the graph and the integral are consistent: We construct a polygon from the graph:

Graphics@Cases[plot, 
  Line[p_] :> Polygon@Reverse@Join[p, {{2 Pi, -2}, {0, -2}}], 
  Infinity]

Then we can take the area of the polygon and subtract the area of the rectangle below the x-axis and approximate the integral. (To get the area the edge of the polygon must trace the boundary counterclockwise; hence the Reverse[].)

Graphics`Mesh`MeshInit[];

Cases[plot, 
  Line[p_] :> PolygonArea@Reverse@Join[p, {{2 Pi, -2}, {0, -2}}], 
  Infinity] - 4 Pi
(*  {1.430715867250397`}   *)

So the approximation agrees with the integral to 14 digits, which is pretty good.

Caveat: Somehow using the built-in Area function on the polygon crashes my kernel (V10.4.1). Hence I used the method from How do I calculate the area of a polygon given its coordinates?.

---

Normalization (response to comment)

@bbgodfrey reports in a comment & his answer a difference between V10.4.1 and V11. The ratio between the integrals have a suspiciously familiar value:

2.023337783/1.43071586725
(*  1.41421  *)

If I check DLMF 28.2.30, I see in V10.4.1, the normalization is satisfied:

NIntegrate[MathieuC[-2 q1 b, q1, x]^2, {x, 0, 2 Pi}, 
 PrecisionGoal -> 8, WorkingPrecision -> 100]
(*  3.14159265358979323846...  *)

NIntegrate[MathieuC[-2 q2 b, q2, x]^2, {x, 0, 2 Pi}, 
 PrecisionGoal -> 8, WorkingPrecision -> 100]
(*  3.141592653589793238462...  *)

The documentation page MathieuC does not indicate an update since V3.0.

Answer 2

It also is possible to solve this problem from first principles, obtaining the symmetric Mathieu functions by solving Mathieu's equation.

s = ParametricNDSolveValue[{c''[x] - 2 q (49/100 + Cos[2 x]) c[x] == 0, c[0] == 1, 
    c'[0] == 0}, c, {x, 0, Pi}, {q}, WorkingPrecision -> 60, MaxSteps -> 50000];

q1 = q /. FindRoot[s[q][Pi] == 1, {q, 606}, WorkingPrecision -> 60]
norm1 = Sqrt[2 NIntegrate[(s[q1][x])^2, {x, 0, Pi}, WorkingPrecision -> 60]/Pi];
Plot[s[q1][x]/norm1, {x, 0, Pi}]
(* 605.789297872353700695280254027186743825904841833249641547632 *)

q2 = q /. FindRoot[s[q][Pi] == 1, {q, 1036}, WorkingPrecision -> 60]
norm2 = Sqrt[2 NIntegrate[(s[q2][x])^2, {x, 0, Pi}, WorkingPrecision -> 60]/Pi];
Plot[s[q2][x]/norm2, {x, 0, Pi}]

(* 1036.13835973661722165073772712373596893094034394911310331033 *)

2 NIntegrate[s[q1][x] s[q2][x], {x, 0, Pi}, WorkingPrecision -> 60]/(norm1 norm2)

(* 1.43127146207803138094684879314159579110898733451406636889128 *)

Results are about the same with WorkingPrecision -> 30.

Although values of q1 and q2 obtained here agree with those in the answer by Michael E2 to of order 10^-29, the final integral differs by of order 0.1% between the two solutions. To attempt to understand this difference, I ran the first block of Michael E2's code on version 10.4.1 and obtained his final integral. However, running his code on version 11.0 produced a final integral of 2.023337783...! The normalization of MathieuC[-2 q1 b, q1, x] differs substantially between the two versions. As noted in the comments attached to the original question, Mathematica's implementation of MathieuC is not reliable.