6
$\begingroup$

So I'm facing the somewhat bizarre situation where I can plot a function and it shows up fine, but trying to numerically integrate it causes problems. This is the plot command.

Plot[MathieuS[MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 1/( 4 0.0275272807115071^2), x] MathieuCPrime[MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x], {x, 0, 6.28}, PlotRange -> Full]

enter image description here

This is what the plot looks like. But when I try to numerically integrate the same function by just replacing 'plot' above with 'NIntegrate' I get this error message:

NIntegrate::inumr: "The integrand MathieuCPrime[-623.771,329.924,x]\ MathieuS[-552.13,329.924,x] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,6.28}}."

Taking the above expression and substituting a value for x doesn't produce a number:

expr = MathieuS[MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x] MathieuCPrime[ MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x]; expr /. x -> 0.1 (* 6.98552*10^-13 MathieuS[-552.13, 329.924, 0.1] *)

Any clues how to proceed?

Improve this question
$\endgroup$
1

5 Answers 5

Reset to default
7
$\begingroup$

This is not a full answer, just a long comment to give more insight into the difference between evaluating the expression within Plot or NIntegrate.

(* note the use of SetDelayed *) expr := MathieuS[ MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x] MathieuCPrime[ MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x]

Substitute in a value for x first, evaluate afterwards:

Block[{x = 0.1}, expr] (* -2.09825*10^-25 *)

Evaluate symbolically first, substitute value afterwards:

Block[{x = 0.1}, Evaluate[expr]] (* 6.98552*10^-13 MathieuS[-552.13, 329.924, 0.1] *)

Not being familiar with this function, I do not know if there is any reason, precision related or other, which makes evaluation difficult for these particular numerical values. If there are precision related reasons, that could be important for whether the result is trustworthy ...

But this observation above does give us a clear workaround for the problem. Just prevent symbolic evaluation within NIntegrate.

f[x_?NumericQ] := MathieuS[MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x] MathieuCPrime[ MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x] NIntegrate[f[x], {x, 0, 6.28}] (* -13.2952 *)
Improve this answer
$\endgroup$
6
$\begingroup$

Use exact numbers in expr by using Rationalize

expr = MathieuS[ MathieuCharacteristicB[2 0 + 2, 1/(4 Rationalize[0.0275272807115071, 0]^2)], 1/(4 Rationalize[0.0275272807115071, 0]^2), x] MathieuCPrime[ MathieuCharacteristicA[0, 1/(4 Rationalize[0.0275272807115071, 0]^2)], 1/(4 Rationalize[0.0275272807115071, 0]^2), x]; Plot[expr, {x, 0, 2 Pi}, PlotRange -> Full]

enter image description here

expr /. x -> 0.1 (* -2.0993*10^-25 *) NIntegrate[expr, {x, 0, 2 Pi}] (* -13.2952 *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ any thoughts on why this is necessary? I assume you've seen this sort of thing before since it seems like one of those tricks everyone tends to accrue. $\endgroup$ Commented Nov 8, 2016 at 16:03
  • 2
    $\begingroup$ @N.J.Evans - Mathematica's internal algorithms generally operate in the complex domain. Consequently, complicated evaluations can accumulate round-off errors that result in complex artifacts. Exact numbers facilitate cancellation of the complex components in what should be real-valued expressions. $\endgroup$ Commented Nov 8, 2016 at 17:45
4
$\begingroup$

I don't know why this happens but a short and ugly hack would be to first evaluate the function on many points of the interval and then integrate the list of evaluated points. This can be automated with FunctionInterpolation:

f[x_] := MathieuS[ MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x] MathieuCPrime[ MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x] fi = FunctionInterpolation[f[x], {x, 0, 6.28}] NIntegrate[fi[x], {x, 0, 6.27}]

Note that FunctionInterpolation shouts out some warning messages too, but it works nonetheless.

** EDIT **

Just to add something, regarding Feyre's answer: Using Riemann sums is of course the straightforward way to go. The problem is that sampling the function at regular intervals may be oversampling it in "non-interesting" regions and undersampling in "interesting" ones. This is why I suggested to use FunctionInterpolation, as it automatically samples more points in more interesting regions.

To see this, here's a plot of the function with the interpolated grid also the density of sampled points (i.e. $1/\Delta x$):

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]; xs = First@InterpolatingFunctionCoordinates[fi]; ys = InterpolatingFunctionValuesOnGrid[fi]; Plot[fi[x], {x, 0, 6.28}, Epilog -> {Red, Point@Transpose@{xs, ys}}, PlotRange -> All] ListPlot[Transpose[{MovingAverage[xs, 2], 1/Differences[xs]}]]

enter image description here

Improve this answer
$\endgroup$
3
$\begingroup$

One way around this is to use Riemann Sums explicitly,

With[{Δx = 1/100}, Total[Table[Δx (MathieuS[ MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x] MathieuCPrime[ MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x]), {x, 0, 2 Pi, Δx}]]]

-13.2952

Improve this answer
$\endgroup$
2
  • $\begingroup$ Note that the difference between this approach and mine is that FunctionInterpolation uses a "smarter" sampling of the functions, supposedly evaluating it on a finer grid at points where the function changes more rapidly. $\endgroup$ Commented Nov 8, 2016 at 15:03
  • $\begingroup$ @yohbs you do have two iterations of numerical approximation. But yes, Interpolation is the superior option in most cases, though I've encountered situations where it would fail. $\endgroup$ Commented Nov 8, 2016 at 15:08
3
$\begingroup$

Also, the problem can be solved simply by adding the pattern _Real in the function definition:

f[x_Real] := MathieuS[MathieuCharacteristicB[2 0 + 2, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x] MathieuCPrime[ MathieuCharacteristicA[0, 1/(4 0.0275272807115071^2)], 1/(4 0.0275272807115071^2), x] NIntegrate[f[x], {x, 0, 6.27}] (*-13.2952*)
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ The robust solution is _?NumericQ, not _Real. None of 1, Pi, Sqrt[2], 1+I are Real, yet it's not 100% clear (though I admit likely) that NIntegrate could never plug these in. $\endgroup$ Commented Nov 8, 2016 at 15:26

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.