1
$\begingroup$

Hey guys I have a pretty simple code that consists of an RMS calculation. The issue is that I am getting quite a few error messages when I run it. The following equations are highly published physics equations with confirmed results, so I am at a loss of why mathematica will not give me an answer:

lambda[z_] := Exp[-(z^2/2^2)]
fo[x_?NumericQ] := 
NIntegrate[
1/(x - xprime)^(1/3) lambda'[xprime], {xprime, -Infinity, x}]
rms[y_] := 
Sqrt[NIntegrate[
lambda[x]*(fo[x])^2, {x, -y, y}] - (NIntegrate[
lambda[x]*(fo[x]), {x, -y, y}])^2]

rms[Infinity]
|improve this question
$\endgroup$
  • $\begingroup$ I get 0. + 1.52231 I as answer, with no errors. What result are you expecting? And what Mathematica version are you using? $\endgroup$ – Dr. belisarius Jun 10 '14 at 18:28
  • $\begingroup$ Sorry my apologies. When rms[0] to rms[100] is when I get the errors $\endgroup$ – user1886681 Jun 10 '14 at 18:29
  • $\begingroup$ !Mathematica graphics $\endgroup$ – Dr. belisarius Jun 10 '14 at 18:42
  • $\begingroup$ Interesting, do you mind explaining the plot a bit? Thank you. $\endgroup$ – user1886681 Jun 10 '14 at 18:46
  • $\begingroup$ On a Mac with v9.0.1 with fo restricted to numeric arguments, fo[x_?NumericQ], then for Table[rms[x], {x, 0, 100, 10}] I get {0., 0. + 1.52231 I, 0. + 1.52231 I, 0. + 1.52231 I, 0. + 1.52231 I, 0. + 1.52231 I, 0. + 1.52231 I, 0. + 1.52231 I, 0. + 1.52231 I, 0. + 1.52231 I, 0. + 1.52231 I} $\endgroup$ – Bob Hanlon Jun 10 '14 at 19:47
3
$\begingroup$

I can reproduce an error when I try to evaluate 

rms[500]

This results in the following error message

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one 
of the following: singularity, value of the integration is 0, highly
oscillatory integrand, or WorkingPrecision too small.`

The problem that you are facing is that, for large values of x the expression

1/(x - xprime)^(1/3) lambda'[xprime]

is strongly peaked around xprime=0 due to the exponential in lambda. In addition there is a singularity at the upper limit. The integral over this expression, on the other hand, is almost zero and the singularity doesn't contribute to the integral. Mathematica has problems numerically integrating such a function because it doesn't know where it should refine the sampling for the integral.

If you look at your functions, you will notice that the Gaussian function lambda limits the interesting range to a couple of standard deviations. The value of rms doesn't change, once the argument goes beyond that range. This means that you can effectively, without noticeably losing accuracy, limit the argument of rms.

In the code below I have limited the integration range to 20σ.

sigma = 2;
lambda[z_] := Exp[-(z^2/sigma^2)]
fo[x_?NumericQ] := 
 NIntegrate[
  1/(x - xprime)^(1/3) lambda'[xprime], {xprime, -Infinity, x}] 
rms[y_] := 
 Sqrt[NIntegrate[
    lambda[x]*(fo[x])^2, {x, Max[-y, -20 sigma], 
     Min[y, 20 sigma]}] - (NIntegrate[
      lambda[x]*(fo[x]), {x, Max[-y, -20 sigma], Min[y, 20 sigma]}])^2]

rms[300]
|improve this answer
$\endgroup$
  • $\begingroup$ This approach is the right direction, but same results can be achieved by forcing NIntegrate to start with a finer grid of regions using the options MinRecursion and MaxRecursion. The explanation in this answer is similar to the one in the section "Examples of Pathological Behavior" of the advanced NIntegrate documentation: reference.wolfram.com/language/tutorial/… . $\endgroup$ – Anton Antonov Sep 24 '15 at 12:11

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.