Tag: NIntegrate

numerical integration – NIntegrate results different when using 0 or 0.0

i was trying to compute a convolution numerically, but got unexpected results. In the end i was able to reduce it to the code below. As you can see, the results are different, depending on whether i write (0-x)^2, x^2 or (0.0-x)^2 (x is a real value). i have Mathematica 10.3.
can somebody explain that to me?
thank you very, very much for your help.
all my best wishes
esra

In(77):= NIntegrate(
 0.16076649110325864` E^(-9774.602659078126` Sqrt(    x^2)) (Erfc(-123.28828005937953` Sqrt(10^-4) + 
      10 Sqrt(110/7) Sqrt(x^2/10^-4)) - 
    E^(19549.20531815625` Sqrt(x^2))
      Erfc(123.28828005937953` Sqrt(10^-4) + 
       10 Sqrt(110/7) Sqrt(x^2/10^-4))) ((
    3 (0 - x)^2)/(0.000025` + (0 - x)^2)^(5/2) - 
    1/(0.000025` + (0 - x)^2)^(3/2)), {x, -Infinity, Infinity})

Out(77)= -410.318

In(84):= NIntegrate(
 0.16076649110325864` E^(-9774.602659078126` Sqrt(    x^2)) (Erfc(-123.28828005937953` Sqrt(10^-4) + 
      10 Sqrt(110/7) Sqrt(x^2/10^-4)) - 
    E^(19549.20531815625` Sqrt(x^2))
      Erfc(123.28828005937953` Sqrt(10^-4) + 
       10 Sqrt(110/7) Sqrt(x^2/10^-4))) ((
    3 (0.0 - x)^2)/(0.000025` + (0.0 - x)^2)^(5/2) - 
    1/(0.000025` + (0.0 - x)^2)^(3/2)), {x, -Infinity, Infinity})

Out(84)= -205.159

In(78):= NIntegrate(
 0.16076649110325864` E^(-9774.602659078126` Sqrt(    x^2)) (Erfc(-123.28828005937953` Sqrt(10^-4) + 
      10 Sqrt(110/7) Sqrt(x^2/10^-4)) - 
    E^(19549.20531815625` Sqrt(x^2))
      Erfc(123.28828005937953` Sqrt(10^-4) + 
       10 Sqrt(110/7) Sqrt(x^2/10^-4))) ((3 x^2)/(0.000025` + x^2)^(
    5/2) - 1/(0.000025` + x^2)^(3/2)), {x, -Infinity, Infinity})

Out(78)= -410.318

equation solving – NIntegrate error along with FindRoot

I have a set of code where I want to find the root (zs) of an equation containing an integral. However, when I evaluated it, it gave me a value that is exactly at the boundary of the given range for FindRoot which is obviously not what I expect. Also, it gave an error message. As I believe the value must be around 0.9zh, it doesn’t have to be exact but it is around that range.

ClearAll("Global`*")
d = 3;
toroot(sqes_?NumericQ, zl_?NumericQ, zh_?NumericQ) := sqes - NIntegrate((1/z^d) Sqrt(1/((1 - (z/zh)^(d + 1)) (1 - (z/zl)^(2 d)))), {z, 0, zl}, WorkingPrecision -> 20)
zs(sqes_?NumericQ, zh_?NumericQ) := zl /. FindRoot(toroot(sqes, zl, zh), {zl, 0.9 zh, 0.5 zh, zh}, WorkingPrecision -> 20)

In(10):= zs(0.00002, 10) // Timing

During evaluation of In(10):= NIntegrate::precw: The precision of the argument function (Sqrt(1/((1-z^4/10000) (1-1.88168*10^-6 z^6)))/z^3) is less than WorkingPrecision (20.`).

During evaluation of In(10):= 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.

During evaluation of In(10):= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {3.652050193603170285401141795380844675863135013439961680436877937000144*10^-224}. NIntegrate obtained 2.531979573217887700978365233494695471156925907058695101296135785961228`70.*^78738 and 2.531979573217887700978365233494695471156925907058695101296135785961228`70.*^78738 for the integral and error estimates.

During evaluation of In(10):= NIntegrate::precw: The precision of the argument function (Sqrt(1/((1-z^4/10000) (1-1.88168*10^-6 z^6)))/z^3) is less than WorkingPrecision (20.`).

During evaluation of In(10):= 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.

During evaluation of In(10):= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {3.652050193603170285401141795380844675863135013439961680436877937000144*10^-224}. NIntegrate obtained 2.531979573217887700978365233494695471156925907058695101296135785961228`70.*^78738 and 2.531979573217887700978365233494695471156925907058695101296135785961228`70.*^78738 for the integral and error estimates.

During evaluation of In(10):= NIntegrate::precw: The precision of the argument function (Sqrt(1/((1-z^4/10000) (1-1.88168*10^-6 z^6)))/z^3) is less than WorkingPrecision (20.`).

During evaluation of In(10):= General::stop: Further output of NIntegrate::precw will be suppressed during this calculation.

During evaluation of In(10):= 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.

During evaluation of In(10):= General::stop: Further output of NIntegrate::slwcon will be suppressed during this calculation.

During evaluation of In(10):= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in z near {z} = {3.652050248022958908304522865892510547317085952963167041802578903669435*10^-224}. NIntegrate obtained 2.531979497759017867516866456547902384642408632229175466570199932032303`70.*^78738 and 2.531979497759017867516866456547902384642408632229175466570199932032303`70.*^78738 for the integral and error estimates.

During evaluation of In(10):= General::stop: Further output of NIntegrate::ncvb will be suppressed during this calculation.

During evaluation of In(10):= FindRoot::reged: The point {10.} is at the edge of the search region {5.0000000000000000000,10.000000000000000000} in coordinate 1 and the computed search direction points outside the region.

Out(10)= {0.28125, 10.}

equation solving – NSolve and NIntegrate, or a better approach

I need to define and plot the following function

$$ a(t) := expleft(int Z(t); dtright) $$

where $Z(t)$ is the solution to the equation

$$ 0 = t – 2 int^Z_1 F(x); dx $$

with $F = F(x)$ being a known (but complicated and non-integrable) function.

How do I define and plot the function $a(t)$ in Mathematica?

Here is my attempt with a particular function $F(x)$ that I need to work with:

A = 0;
F(x_) = - ((4*A*x^(9/2) + 64*x^6 - 4*Sqrt(A*x^9*(32*x^(3/2) + A)))^(1/3)/(16*x^4 - 4*x^2*(4*A*x^(9/2) + 64*x^6 - 4*Sqrt(A*x^9*(32*x^(3/2) + A)))^(1/3) - (4*A*x^(9/2) + 64*x^6 - 4*Sqrt(A*x^9*(32*x^(3/2) + A)))^(2/3)));

A = 0;
Int(Z_?NumericQ) := NIntegrate(F(x), {x, 1, Z})
S(t_?NumericQ) := NSolve(t - 2*Int(Z) == 0, Z)
a(t_) := Exp(Integrate(S(t), t))

However, when trying to evaluate for example $a(2)$ I get the following error:

The error

numerical integration – Why does the Kernel CRASH without any Message when I use “NIntegrate” for the following complicated compiled Integrand?

When I NIntegrate the Integrand with the particular parameter showing below the Kernel crash without any Message.

But if I change the “range” parameter a little bit, the Kernel won’t crash. The crash problem seems to occur only under certain parameters.

Because I want to make a density plot of the integral as a function of the parameters, this crash problem under certain parameters always breaks the density plot process.

PS: This crash problem happens both on Mac and Win10

(*Integrand*)

expression1 = {{-0.03783737745791316` + 
     52.07594879537986` (k1^2 + k2^2) - 
     0.0006566020385220937` (GothicCapitalB)3, 
-0.00020215783672119025` ((GothicCapitalB)1 - 
       I (GothicCapitalB)2), (0.` + 0.` I) + 
     1/2 (273.35324964239567` (k1 - I k2)^3 - 
        22.668570452264646` (k1 + I k2)^3), 
    0.` + I (k1 - I k2) (-1.0301133331497063` - 
        370.16816003256605` (k1^2 + 
           k2^2))}, {-0.00020215783672119025` ((GothicCapitalB)1 + 
       I (GothicCapitalB)2), -0.03783737745791316` + 
     52.07594879537986` (k1^2 + k2^2) + 
     0.0006566020385220937` (GothicCapitalB)3, 
    0.` - I (k1 + I k2) (-1.0301133331497063` - 
        370.16816003256605` (k1^2 + k2^2)), (0.` + 0.` I) + 
     1/2 (22.668570452264646` (k1 - I k2)^3 - 
        273.35324964239567` (k1 + I k2)^3)}, {(0.` + 0.` I) + 
     1/2 (-22.668570452264646` (k1 - I k2)^3 + 
        273.35324964239567` (k1 + I k2)^3), 
    0.` + I (k1 - I k2) (-1.0301133331497063` - 
        370.16816003256605` (k1^2 + k2^2)), -0.03304414880767263` + 
     95.54303239292248` (k1^2 + k2^2) - 
     0.0015031216532114644` (GothicCapitalB)3, 
-0.00017807130082869512` ((GothicCapitalB)1 - 
       I (GothicCapitalB)2)}, {0.` - 
     I (k1 + I k2) (-1.0301133331497063` - 
        370.16816003256605` (k1^2 + k2^2)), (0.` + 0.` I) + 
     1/2 (-273.35324964239567` (k1 - I k2)^3 + 
        22.668570452264646` (k1 + 
           I k2)^3), -0.00017807130082869512` ((GothicCapitalB)1 + 
       I (GothicCapitalB)2), -0.03304414880767263` + 
     95.54303239292248` (k1^2 + k2^2) + 
     0.0015031216532114644` (GothicCapitalB)3}};
expression2 = {{104.15189759075972` k1s, 0, 
    1/2 (820.059748927187` (k1s - I k2)^2 - 
       68.00571135679394` (k1s + I k2)^2), (0.` - 
        740.3363200651321` I) k1s (k1s - I k2) + 
     I (-1.0301133331497063` - 
        370.16816003256605` (k1s^2 + k2^2))}, {0, 
    104.15189759075972` k1s, (0.` + 740.3363200651321` I) k1s (k1s + 
        I k2) - I (-1.0301133331497063` - 
        370.16816003256605` (k1s^2 + k2^2)), 
    1/2 (68.00571135679394` (k1s - I k2)^2 - 
       820.059748927187` (k1s + I k2)^2)}, {1/
     2 (-68.00571135679394` (k1s - I k2)^2 + 
       820.059748927187` (k1s + I k2)^2), (0.` - 
        740.3363200651321` I) k1s (k1s - I k2) + 
     I (-1.0301133331497063` - 370.16816003256605` (k1s^2 + k2^2)), 
    191.08606478584497` k1s, 
    0}, {(0.` + 740.3363200651321` I) k1s (k1s + I k2) - 
     I (-1.0301133331497063` - 370.16816003256605` (k1s^2 + k2^2)), 
    1/2 (-820.059748927187` (k1s - I k2)^2 + 
       68.00571135679394` (k1s + I k2)^2), 0, 
    191.08606478584497` k1s}};
expression3 = {{104.15189759075972` k2s, 0, 
    1/2 ((0.` - 820.059748927187` I) (k1 - I k2s)^2 - (0.` + 
          68.00571135679394` I) (k1 + 
          I k2s)^2), -1.0301133331497063` - (0.` + 
        740.3363200651321` I) (k1 - I k2s) k2s - 
     370.16816003256605` (k1^2 + k2s^2)}, {0, 
    104.15189759075972` k2s, -1.0301133331497063` + (0.` + 
        740.3363200651321` I) (k1 + I k2s) k2s - 
     370.16816003256605` (k1^2 + k2s^2), 
    1/2 ((0.` - 68.00571135679394` I) (k1 - I k2s)^2 - (0.` + 
          820.059748927187` I) (k1 + I k2s)^2)}, {1/
     2 ((0.` + 68.00571135679394` I) (k1 - I k2s)^2 + (0.` + 
          820.059748927187` I) (k1 + 
          I k2s)^2), -1.0301133331497063` - (0.` + 
        740.3363200651321` I) (k1 - I k2s) k2s - 
     370.16816003256605` (k1^2 + k2s^2), 191.08606478584497` k2s, 
    0}, {-1.0301133331497063` + (0.` + 740.3363200651321` I) (k1 + 
        I k2s) k2s - 370.16816003256605` (k1^2 + k2s^2), 
    1/2 ((0.` + 820.059748927187` I) (k1 - I k2s)^2 + (0.` + 
          68.00571135679394` I) (k1 + I k2s)^2), 0, 
    191.08606478584497` k2s}};
HamiltCompiled = 
  Compile({{k1, _Real}, {k2, _Real}, {(GothicCapitalB)1, _Real}, {
(GothicCapitalB)2, _Real}, {(GothicCapitalB)3, _Real}}, 
   Evaluate(expression1), 
   CompilationOptions -> {"ExpressionOptimization" -> True});
Dk1HamiltCompiled = 
  Compile({{k1s, _Real}, {k2, _Real}, {(GothicCapitalB)1, _Real}, {
(GothicCapitalB)2, _Real}, {(GothicCapitalB)3, _Real}}, 
   Evaluate(expression2), 
   CompilationOptions -> {"ExpressionOptimization" -> True});
Dk2HamiltCompiled = 
  Compile({{k1, _Real}, {k2s, _Real}, {(GothicCapitalB)1, _Real}, {
(GothicCapitalB)2, _Real}, {(GothicCapitalB)3, _Real}}, 
   Evaluate(expression3), 
   CompilationOptions -> {"ExpressionOptimization" -> True});
OccupiedCurvatureCompiled((Mu)_?NumericQ, k1_?NumericQ, 
   k2_?NumericQ, (GothicCapitalB)1_, (GothicCapitalB)2_, 
(GothicCapitalB)3_) := 
  Block({Mx1, Mx2, values, vectors, curvaturesList},
   {values, vectors} = 
    Eigensystem@
     HamiltCompiled(k1, 
      k2, (GothicCapitalB)1, (GothicCapitalB)2, 
(GothicCapitalB)3);
   {values, vectors} = {values((#)), vectors((#))} &(
     Ordering@values);
   Mx1 = vectors(Conjugate).(Dk1HamiltCompiled(k1, 
       k2, (GothicCapitalB)1, (GothicCapitalB)2, 
(GothicCapitalB)3)).vectors(Transpose);
   Mx2 = vectors(Conjugate).(Dk2HamiltCompiled(k1, 
       k2, (GothicCapitalB)1, (GothicCapitalB)2, 
(GothicCapitalB)3)).vectors(Transpose);
   curvaturesList = 
    Re(I Total((# - #(Transpose)))) &(
     Mx1*Mx2(Transpose)*
      Table(If(n1 == n2, 0, 1/(values((n1)) - values((n2)))^2), {n1, 
        1, 4}, {n2, 1, 4}));
   curvaturesList.Boole@Thread(values < (Mu))
   );

(*No Crash*)
With({range = 
   0.0006` + 
    0.00025*19.`, (Mu) = -0.034637500000000015`, (GothicCapitalB)1 =
    0., (GothicCapitalB)2 = 19.`, (GothicCapitalB)3 = 0.}, 
 NIntegrate(
  OccupiedCurvatureCompiled((Mu), k1, 
   k2, (GothicCapitalB)1, (GothicCapitalB)2, (GothicCapitalB)3), 
{k1, -range, range}, {k2, -range, range}, MinRecursion -> 2, 
  PrecisionGoal -> 3, AccuracyGoal -> 3))

(*Crash*)
With({range = 
   0.0005` + 
    0.00025*19.`, (Mu) = -0.034637500000000015`, (GothicCapitalB)1 =
    0., (GothicCapitalB)2 = 19.`, (GothicCapitalB)3 = 0.}, 
 NIntegrate(
  OccupiedCurvatureCompiled((Mu), k1, 
   k2, (GothicCapitalB)1, (GothicCapitalB)2, (GothicCapitalB)3), 
{k1, -range, range}, {k2, -range, range}, MinRecursion -> 2, 
  PrecisionGoal -> 3, AccuracyGoal -> 3))

numerical integration – ParallelTable of NIntegrate of an Interpolating Function

I construct a 2D interpolating function using the following data (MWE):

Rdata = Flatten(Table({{theta, phi}, 5.*Cos(theta + phi)}, {theta, 0, 2. Pi, 0.1}, {phi, 0, 2. Pi/6, 0.1}), 1);
Rinterp = Interpolation(Rdata, InterpolationOrder -> 1);

Now, since I want to use this interpolating function to do an integral, I use the following to make things faster:

RinterpQ(theta_?NumericQ, phi_?NumericQ) := Rinterp(theta, phi);

And now I use ParallelTable to evaluate the following integrals:

ParallelTable(
 NIntegrate(
  RinterpQ(theta, phi)*Cos(n*5*phi), {theta, 0, 2. Pi}, {phi, 0, 
   2. Pi/5}, Method -> {Automatic, "SymbolicProcessing" -> 0}), {n, 0,
   1})

which in my computer with 4 threads is done in 26s.

However, if I increase a lot the number of data I use for the interpolating function, say

fdata = Flatten(Table({{theta, phi}, 5.*Cos(theta + phi)}, {theta, 0, 2. Pi, 0.001}, {phi, 0, 2. Pi/6, 0.001}), 1);

then, by repeating the interpolating process and defining RinterpQ, the ParallelTable evaluation never ends!

Notice also that if I don’t do the RinterpQ, and I directly use Rinterp directly, I don’t have any problems with the fine grid. And also if I just use Parallel I don’t have any type of problems with the fine grid.

PS: I know that using ParallelTable to make only two evaluations is not smart, but this is just an example!

numerical integration – Help with NIntegrate settings to evaluate integral

I am trying to evaluate this integral:
begin{align*}
alpha_{2}=int_{-infty}^{1.645}left(1-Phileft(frac{sqrt{25}}{sqrt{15}} 1.645-frac{sqrt{10}}{sqrt{15}} z_{1}right)right) phileft(z_{1}right) d z_{1}
end{align*}
where $Phi$ is the Normal CDF and $phi$ is the normal PDF. I know that the answer should be 0.03325.

I used the following code, but it doesn’t converge to an answer. Any suggestions?

pdf(x_) := PDF(NormalDistribution(0, 1), x)

cdf(x_) := CDF(NormalDistribution(0, 1), x)

NIntegrate( 1 - cdf(Sqrt(25)/Sqrt(15) 1.645 - Sqrt(10)/Sqrt(15) x)* pdf(x), {x, -Infinity, 1.645})

which returns

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.

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {-8.16907*10^224}. NIntegrate obtained 8.582639123060543`15.954589770191005*^27949 and 8.582639123060543`15.954589770191005*^27949 for the integral and error estimates.
```

NIntegrate of a convergent integral working with large integration limits, but not with infinite integration limits

Code for reproducing is below:

Integrand(x_) := ((x + 1) Abs(x) Log((x + 1)^2/Abs(x)^3))/((x + 1)^2 - Abs(x)^3);
    
NIntegrate(Integrand(x), {x, -10000000, 10000000})
NIntegrate(Integrand(x), {x, -Infinity, Infinity})

The first instance of NIntegrate gives $2.48398$, and one can check by plotting the result of NIntegrate as a function of the integration limits that the integral appears to converge to this answer for very large integration limits. However, the second instance of NIntegrate gives a completely different answer of $1.75434 times 10^{8}$. What’s going on here, and how can I make the two integrals agree? I think something that NIntegrate might be having trouble with is that the convergence of the integral requires the fact that the integrand is essentially odd for very large $|x|$ and these contributions tend to cancel out. Ideally I’d like to find some way to get Mathematica to deal with this properly with the infinite integration range, so that I don’t have to keep remembering to exchange the infinite integration limits with large finite values every time I need to do an integral like this.

NIntegrate fails to give non-numerical value error, instead returns zero

When I evaluate

NIntegrate(a(x - y)^2, {x, 0, 1}, {y, 0, 1}, AccuracyGoal -> 8)

Mathematica returns 0 without giving a non-numerical value error. If I get rid of the AccuracyGoal specification then the expected error is raised. It is also raised if I multiply a by pretty much any function other than (x-y)^2.

I’d like to understand why this is happening. In the past, I was occasionally sloppy and would append an expression like the one above with /. {a -> 1}; this would give an error but nevertheless return the correct answer. In the example above, doing so just gives the wrong answer without any error.

numerical integration – Shouldn’t NIntegrate return a number whose precision is PrecisionGoal, not WorkingPrecision?

I know that Mathematica has great built-in precision tracking, so when you do calculations with arbitrary-precision numbers, Mathematica keeps track of the precision on the result. Given this careful attention to numerical error and precision tracking, I am surprised that, say

InputForm(NIntegrate(E^(-x^2), {x, 0, Infinity},PrecisionGoal -> 20, WorkingPrecision -> 100))

returns a number with precision 100, not 20. I know Mathematica is using precision-100 numbers in its numerical calculations for NIntegrate, but the function is built to return a number whose actual precision is at least 20. In the spirit of useful precision tracking, wouldn’t it make more sense for NIntegrate to return a number with a precision of PrecisionGoal, not WorkingPrecision?

This question is more about numerical coding philosophy than about how NIntegrate works. But this is important as Wolfram presumably makes these decisions with use cases in mind, so I want to know if I’m missing something.

“Further output of NIntegrate” ERROR at diffusion equation

Good day everybody,

I am trying to solve the diferential equation for diffusion:

dt(u)=A*dx(u^(m)*dx(u)),

where “m” is the difussion coefficient. In this case, I am trying to make “m” a function of x, and solve the equation numerically.

GCo = First(u /. NDSolve({D(u(x,t),t)==dCo*D((u(x,t)^mCo(x)*D(u(x,t),x)),x),(D(u(x,t),x)./->-20)=0,(D(u(x,t),x)./->40)=0,u(x,0)=FC0(x)},u,{x,-20,40},{t,0,10}))

Where “FCo(x)” is the initial conditions.

I get the following errors:

“General::stop: Further output of NIntegrate::nlim will be suppressed during this calculation.”

“General::stop: Further output of General::munfl will be suppressed during this calculation.”

Is there any way to solve this?

Thank you very much for your help