6
$\begingroup$

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!

$\endgroup$
  • 1
    $\begingroup$ I believe this is the cause of the problem: MathieuC[-2 q b, q, x] doesn't seem to evaluate at any x. $\endgroup$ – JungHwan Min Aug 28 '16 at 15:53
  • $\begingroup$ Trouble with Mathieu functions: stackoverflow.com/q/7798435, (58297), (81395), (83032). $\endgroup$ – Michael E2 Aug 28 '16 at 17:31
  • $\begingroup$ @MichaelE2 All the linked issues still persist in v11 on my machine. Can you confirm? $\endgroup$ – QuantumDot Aug 28 '16 at 19:42
  • $\begingroup$ @QuantumDot No, I'm using 10.4.1. $\endgroup$ – Michael E2 Aug 28 '16 at 19:47
  • $\begingroup$ @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. $\endgroup$ – bbgodfrey Aug 29 '16 at 5:30
5
$\begingroup$

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}]

Mathematica graphics

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

Mathematica graphics

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]

Mathematica graphics

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – confused Aug 29 '16 at 10:59
  • $\begingroup$ @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.99`100 in my code above, the computations go through and seem reasonable (i.e. consistent). I'm using V10.4.1. $\endgroup$ – Michael E2 Aug 29 '16 at 12:21
  • $\begingroup$ 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. $\endgroup$ – confused Aug 29 '16 at 13:30
  • $\begingroup$ @Assaf Yes, I got negative -- didn't think about that. Oh well...Good luck! $\endgroup$ – Michael E2 Aug 29 '16 at 16:48
4
$\begingroup$

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 *)

enter image description here

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 *)

enter image description here

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.

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

Your Answer

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

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