Numerical error in Mathieu functions - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-10-19T19:06:54Z https://mathematica.stackexchange.com/feeds/question/58297 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/58297 6 Numerical error in Mathieu functions xslittlegrass https://mathematica.stackexchange.com/users/1364 2014-08-27T17:01:53Z 2015-05-16T17:08:07Z <p>Consider the <a href="http://reference.wolfram.com/mathematica/ref/MathieuCharacteristicA.html" rel="nofollow noreferrer"><code>MathieuCharacteristicA</code></a> function, which is a piecewise function according to the documentation. The discontinuity happens at integer number.</p> <pre><code>With[{V0 = -1}, Plot[MathieuCharacteristicA[κ, V0], {κ, -2.5, 2.5}]] </code></pre> <p><img src="https://i.stack.imgur.com/yn3nY.png" alt="enter image description here"></p> <p>Consider the point approaching k=2 from the left side , and plot the Mathieu funtions near that point. </p> <pre><code>ParallelTable[ Plot[Evaluate@ With[{V0 = -1, κ = 2 - ϵ}, Re@MathieuC[MathieuCharacteristicA[κ, V0], V0, z]], {z, -10, 10}, PlotRange -&gt; All, ImageSize -&gt; Medium], {ϵ, {10^-8, 15/10*10^-8, 18/10*10^-8, 2*10^-8}}] </code></pre> <p><img src="https://i.stack.imgur.com/ZXDAl.png" alt="enter image description here"></p> <p>We see that from points k=2-10^-8 , k=2-1.5*10^-8 to points k=1.8*^-8, k=2*^-8, there are big discontinouity. <strong>Why does this big discontinuity happen in the Mathieu function</strong>, even we are still away from the piecewise point? <strong>Which result is correct?</strong></p> <p>Moreover, as I increase the working precision, the results changes. Which results should I trust?</p> <pre><code>ParallelTable[ Plot[Evaluate@ With[{V0 = -1, κ = 2 - ϵ}, Re@MathieuC[MathieuCharacteristicA[κ, V0], V0, z]], {z, -10, 10}, PlotRange -&gt; All, ImageSize -&gt; Medium, WorkingPrecision -&gt; 50], {ϵ, {10^-8, 15/10*10^-8, 18/10*10^-8, 2*10^-8}}] </code></pre> <p><img src="https://i.stack.imgur.com/YYeZ3.png" alt="enter image description here"></p> <p><strong>Update:</strong></p> <p>More strange behavior</p> <pre><code>NLimit[ Re@MathieuC[MathieuCharacteristicA[κ, -1], -1, 0], κ -&gt; 2, Direction -&gt; 1, WorkingPrecision -&gt; 100] (* 0.000026560352729499428275267693547091828644960849846890155742135607985075453865741662994877041 *) N[ Table[Re@MathieuC[MathieuCharacteristicA[2 - ϵ, -1], -1, 0], {ϵ, {10^-6, 10^-8, 10^-10, 10^-20, 10^-40, 10^-60, 10^-100}}], 100] (* \ {9.375519741470728355592491183508603286638427801561870220416306315833951776806837902179623867179198570*10^-6, 9.375519742493871990285719573456924995106820123921565403331967009687923231481580841077469030562191305*10^-8, 9.375519742493974304649207591451232553957148872416907065490732452953108780592926489536671069329119067*10^-10, 9.375519742493974314881667185336397875216528878100913700967589255793432748631135267810942389900367393*10^-20, 9.375519654864253585910474819580416042344587293945440367769556457939085038181962308922619162634748808*10^-40, 1.114388591781733115021002428520876171768008184684161143978162399223107768160400228809114697590700165, 1.114388591781733115021002428520876171768008184684161143978162399223107768160400228809114697590700165} *) </code></pre> <p>What's the correct limit for k->2 ?</p> https://mathematica.stackexchange.com/questions/58297/-/58301#58301 5 Answer by Michael E2 for Numerical error in Mathieu functions Michael E2 https://mathematica.stackexchange.com/users/4999 2014-08-27T17:35:24Z 2014-08-28T01:14:57Z <p>If we compare <code>N</code> on the exact values with the <code>MachinePrecision</code> values, we see that the second two graphs (of the first quartet) look correct and the first two are wrong.</p> <pre><code>Block[{z = 0}, Table[With[{V0 = -1, κ = 2 - ϵ}, Re@MathieuC[MathieuCharacteristicA[κ, V0], V0, z]], {ϵ, {15/10*10^-8, 18/10*10^-8}}] ] N[%, 6] (* {MathieuC[MathieuCharacteristicA[399999997/200000000, -1], -1, 0], MathieuC[MathieuCharacteristicA[999999991/500000000, -1], -1, 0]} {1.4063279613740616126811177594156`6.*^-7, 1.6875935536488557009743434482017`6.*^-7} *) Block[{z = 0.}, Table[With[{V0 = -1, κ = 2 - ϵ}, Re@MathieuC[MathieuCharacteristicA[κ, V0], V0, z]], {ϵ, {15/10*10^-8, 18/10*10^-8}}] ] (* {1.11439, 1.7027*10^-7} *) </code></pre> <p>As the OP showed, this can be seen in the plots, if the <code>WorkingPrecision</code> is set high enough. Clearly, one should trust the second quartet of plots if one is going to trust <em>Mathematica</em> at all. <code>N[expr, n]</code> will report the answer with a precision that is supposed to be correct. </p> <hr> <p><em>Edit</em></p> <p>As <a href="https://mathematica.stackexchange.com/users/16/acl">@acl</a> has observed Mathieu functions are difficult functions numerically. I would expect it to be even more difficult near a singular point of <code>MathieuCharacteristicA[κ, -1]</code>, where a little round-off error causes a discontinuous jump.</p> <p>The following seem entirely consistent with the left-hand limit being zero, which is at the same time not inconsistent with the OP's evaluation of <code>NLimit</code>.</p> <pre><code>N@Block[{\$MaxExtraPrecision = 500}, NLimit[Re@MathieuC[MathieuCharacteristicA[κ, -1], -1, 0], κ -&gt; 2, Direction -&gt; 1, WorkingPrecision -&gt; 300, Terms -&gt; 50] ] N[Re@MathieuC[MathieuCharacteristicA[2, -1], -1, 0], 10] N@Block[{\$MaxExtraPrecision = 500}, NLimit[Re@MathieuC[MathieuCharacteristicA[κ, -1], -1, 0], κ -&gt; 2, Direction -&gt; -1, WorkingPrecision -&gt; 300, Terms -&gt; 50] ] (* -1.49047*10^-45 0.8157268391 1.15361 *) </code></pre> <p>The inconsistencies observed by the OP seem to be due to round-off error. I suppose one complaint is that <em>Mathematica</em> issues no warnings in evaluating the OP's examples, especially the ones involving <code>N</code> applied to an exact value.</p>