Mathematica doesn't evaluate MathieuC function for large arguments

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!

• I believe this is the cause of the problem: MathieuC[-2 q b, q, x] doesn't seem to evaluate at any x. Aug 28 '16 at 15:53
• Trouble with Mathieu functions: stackoverflow.com/q/7798435, (58297), (81395), (83032). Aug 28 '16 at 17:31
• @MichaelE2 All the linked issues still persist in v11 on my machine. Can you confirm? Aug 28 '16 at 19:42
• @QuantumDot No, I'm using 10.4.1. Aug 28 '16 at 19:47
• @MichaelE2 Running the first block of your code with version 11.0, I do not obtain the same results that you (and I) obtain with version 10.4.1. Aug 29 '16 at 5:30

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.49100;
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[].)

GraphicsMeshMeshInit[];

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.

• Thanks! This does seem to work for b=0.49 as I originally defined. The problem is that in my project I also need to work with values like b=0.99, for which the precision just doesn't seem to be enough (I tried up to WorkingPrecision=10,000) or maybe I am getting something wrong... Do you have an idea how to deal with that? Aug 29 '16 at 10:59
• @Assaf Maybe you can use the diff. eq. like Ruslan did in his question here and bbgodfrey did his answer to this question? This shows how sensitive the built-in MathieuC is to slight errors: Table[MathieuC[-2 q2 b, q, 11/10.100], {q, q2 - 10^-20, q2 + 10^-20, 10^-20}] (use with the 100-digit parameters in my answer). Actually, though, if I change b to b = 0.99100 in my code above, the computations go through and seem reasonable (i.e. consistent). I'm using V10.4.1. Aug 29 '16 at 12:21
• I will try using the diff eq. I am guessing when you changed to b=0.99 you got negative values for q. for them it works, but there are also positive values (they are very large, so the start point for FindRoot should be about 10^7 to find them), and when I try it with them it fails - I get non-periodic functions that grow exponentially. I should have mentioned that earlier. Aug 29 '16 at 13:30
• @Assaf Yes, I got negative -- didn't think about that. Oh well...Good luck! Aug 29 '16 at 16:48

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.

• I looked into this and wondered if MathieuC was discontinuous at the value of q considered. Aug 29 '16 at 5:42