Function providing input and integration limits for NIntegrate

I am trying to define some custom function that evaluates the numeric integral of some complicated function. The problem is that I would like to include as an input the integration limits. A MWE follows

f[x_,y_]:= Exp[-2 x^2 y^2]
test[x_?NumericQ,IntL_]:= NIntegrate[f[x, y],IntL]
myIntL1= {y,-4 x^2,4x^2};
myIntL2= {y,-4 x^4,4x^4};

Then, if I evaluate for instance test[3,myIntL1] I get a problem concerning an invalid limit of integration.

Is there a clever way to fix this without defining several functions including the integration limits, such as

    test1[x_?NumericQ,IntL_]:= NIntegrate[f[x, y],{y,-4 x^2,4x^2}]
    test2[x_?NumericQ,IntL_]:= NIntegrate[f[x, y],{y,-4 x^4,4x^4}]


Of course, here everything looks simple, but in my case all the functions are rather long; some even purely numeric ones. Since I have multiple choices for the integration limits, it would be more practical to avoid defining test1[x_], test2[x_], ..., testN[x_]

Thanks in advance,

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 – Help with NIntegrate settings to evaluate integral

I am trying to evaluate this integral:
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}

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.

The following code in R gives me the correct answer:

inside <- function(z1, n1, n2, cv) {
  nt <- n1 + n2
  (1 - pnorm(sqrt(nt/n2) * cv - sqrt(n1/n2) * z1)) * dnorm(z1)

additional.error <- function(n1, n2, cv) {
  integrate(inside, lower = -Inf, upper = cv, n1 = n1, n2 = n2, cv = 1.645)$value

additional.error(n1 = 10, n2 = 15, cv = qnorm(0.95))

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:


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

differential equations – Solve ODE via NDSolveValue with NIntegrate of the Undetermined Function Being one of Terms

I’m new here. If there is anything not appropriate pls let me know.

I am currently working on a differential equation with one of which term is a integral of the variable.

frac{d^2u(x)}{dx^2}=cosh(G(x))+frac{1}{C_1}int_{0}^{1}{u(x)sinh(G(x))dx }+C_2

with the boundary condition, $ u(x=0)=0$ and $u'(x=1)=0$,

where $ G(x)= 2ln(frac{1+C_3cdot exp(-x)}{1-C_3cdot exp(-x)})$ and
C1, C2 and C3 are system constants which could be predefined.

To show it more simply, I let $C_1=C_2=C_3=1$ in the code below,

G(x_) = 2 Log((1 + Exp(-x))/(1 - Exp(-x)));
Sol = 
    u''(x) == 
       Cosh( G(x) ) + NIntegrate( u(x) *Sinh( G(x) ),{x,0,1}) + 1
       ,u'(1) == 0., u(0) == 0.
  , u, {x, 0, 1}, PrecisionGoal -> 10) ;

However, since u(x) in the integral have yet been solved, the numerical integration would fail with the error message show:

"The integrand {} has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1}}"*

Some suggested by breaking the procedure of NDSolveValue in parts with the NIntegrate inserted. However, I am not sure how to do it correctly in Mathematica.

Thanks for your kindly help, I am really appreciated it!

numerical integration – How to solve the problem in NIntegrate NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy

The integration is actually NIntegrate(-0.17116940940118283` + 1/( 9.736942322213501` + 7.789553857770802` Cos(q)) + ( 0.02866566930866079` (0.5` + 1.` Cos(q)) Sin( q) (-3.0095696738628313` Sqrt(1.25` + 1.` Cos(q)) Cos(0.` + ArcTan((0.5` Sin(q))/(-1 - 0.5` Cos(q)))) + 1.` Sin(q)))/( 0.9772727272727273` + 1.` Cos(q) - 0.045454545454545456` Cos(2 q) - 0.09090909090909091` Cos(3 q)) + ((0.35586923225834494` + 0.5931153870972414` Cos(q) + 0.11862307741944829` Cos(2 q)) Sin( 0.` + ArcTan((0.5` Sin(q))/(-1 - 0.5` Cos(q)))))/((1.75` + 1.` Cos(q) - 0.5` Cos(2 q))^(3/2) Sqrt( 1 - (1.` Sin(q)^2)/( 1.75` + 1.` Cos(q) - 0.5000000000000001` Cos(2 q)))), {q, -Pi, Pi}). Error message is **NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in q near {q} = {-3.14159254089972008785892145083358745552559732061581598827615380287}. NIntegrate obtained -1.24910^-16 and 4.588053980254483`^-13 for the integral and error estimates.**How to get the real integration value?

numerical integration – Problem with NIntegrate over a piecewise function

The following code seems to generate wrong results

q2 = 1/y^0.9;
G2 = 1.1 x^0.045;
b2bar = x /. Solve[G2 == 1, x][[1]]
ff = FullSimplify[PiecewiseExpand[D[Min[G2, 1], x]*x*q2]]
NIntegrate[ff, {y, 0, 1}, {x, y, b2bar}]
NIntegrate[ff, {y, 0, 1}, {x, y, 1}]

Since ff=0 for x>=b2bar, one would expect NIntegrate[ff, {y, 0, 1}, {x, y, b2bar}] and NIntegrate[ff, {y, 0, 1}, {x, y, 1}] to generate the same results. But NIntegrate[ff, {y, 0, 1}, {x, y, b2bar}] gives 0.038246 and NIntegrate[ff, {y, 0, 1}, {x, y, 1}] gives 0.022375.

Is there an error with NIntegrate when integrating piecewise functions? How can I avoid this?

numerical integration – NIntegrate blowing up/behaving weirdly at the turning point of the integrand

I’m performing (what should be) a straightforward numerical integration (Fourier transform) of a function with no poles / singularities (at least in a particular parameter regime):

<< FourierSeries`


l = 1;
ϵ = 10^-4;
r = 0.1; (*radial coordinate*)

p = 0;
P = 80;

Data1 = ConstantArray(Null, {P, 2});

While(p < P, p++;
  w = -10 + 20 p/P;
  σ = (
   2 (-2 l^2 r^2 Sin(s/(2 Sqrt(l^-2 - r^2)))^2 + (1 - l^2 r^2) 2 Sinh(
        s/(2 Sqrt(l^-2 - r^2)) - I ϵ)^2))/l^2;
  RF = NInverseFourierTransform(-1/(4 π^2) E^(-I w s)/σ, s,
      w, Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
     MinRecursion -> 6) // Chop;
  Data1((p, 1)) = w;
  Data1((p, 2)) = RF;

ListLinePlot(Data1, PlotRange -> All, 
   LabelStyle -> {FontSize -> 16, FontFamily -> "CMU Serif", Black}, 
   PlotLabel -> StringForm("r = ``.", {r // N})) // Print;

When the parameter r is small, the numerics seem to work fine/as expected, producing plots like this (r = 0.1):
enter image description here

(x-axis is the energy, y-axis is the transition rate of a two-level system – this shows a roughly thermal Planck spectrum).

When increasing r beyond some critical point, the numerics seem to blow up, giving a weird answer. For example, r = 0.8 yields:

enter image description here

(see especially the magnitude of the y-axis).

I plotted the integrand for these two values below. The numerics seem to blow up when the function bifurcates from having one turning point to two turning points:

r=0.1 (integrand plotted as a function of s)

enter image description here

and r =0.8

enter image description here

Why does NIntegrate not like this seemingly inconspicuous change in behaviour of the integrand? Any help would be greatly appreciated! (If there’s an analytic solution, even better :P)

equation solving – Nested NIntegrate with FindRoot

I am trying to numerically integrate a function using nested NIntegrate:

$$F(N,x,s)=int_{-infty}^s int_{-infty}^{+infty} K(N,z’,x,x’) g_{x’,s’} dx’ds’ $$

where the kernel of the integration, $K(N,z’x,x’)$, is a messy expression defined in the mathematica code below, and $g_{x’,s’}$ is a bi-variate gaussian defined by:

$$ g_{x,s’}=frac{n}{2pisigma_{x’}sigma_{s’}}expleft({ -frac{x’^2}{2sigma_{x’}^2} }right)expleft({ -frac{s’^2}{2sigma_{s’}^2} }right).$$

The tricky part(s) is that:

  1. $z’$ in the $K(N,z’,x,x’)$ needs to be solved for numerically using FindRoot and will have a $s’$ dependence.
  2. The integration upper limit over $ds’$ is a variable $s$.
  3. I suspect the kernel is oscillatory with $N$ (denoted “Kernel” in the code below) so maybe an averaging of the kernel over $N$ can be done to simplify the kernel and eliminate $N$ if the integrations prove to be too time consuming.

At the end, I would like a function, F(N,x,s), that would be able to plot across $s$ for a given $(N,x)$ values i.e. Plot(F(a,b,s,{s,-1e-5,1e-5}).

e = -1.60217733*10^-19;
m = 9.109389699999999*10^-31;
epsilon = 8.854187817620391*10^-12;
re = 2.81794092*10^-15;
c = 2.99792458*10^8;
n = -10^-10/e;
KK = 1;
lw = 0.026;
kw = (2 Pi)/lw;
gamma = 4000/0.511;
beta = Sqrt(1 - 1/gamma^2);
sigmaS = 10^-5;
sigmaX = 30*10^-6;
coeff = n/(2 Pi*sigmaS*sigmaX) Exp(-(xprime^2/(2 sigmaX^2)))*
Exp(-(sprime^2/(2 sigmaS^2)));

(*Preliminary Equations*)
rs2 = {zprime, xprime + KK/(gamma*kw) Sin(kw*zprime), 0};
ro2 = {(NN + 10000)*lw, x + KK/(gamma*kw) Sin(kw*(NN + 10000)*lw), 0};

betas = {beta - KK^2/(2 gamma^2) Cos(kw*zprime)^2,KK/gamma Sin(kw*zprime), 0};
betao = {beta - KK^2/(2 gamma^2) Cos(kw*(NN + 10000)*lw)^2,KK/gamma Sin(kw*(NN + 10000)*lw), 0};

betaDot = {(c*KK^2*kw)/(2 gamma^2)Sin(2 kw*zprime), -((KK*c*kw)/gamma) Sin(kw*zprime), 0};

deltar2 = ro2 - rs2;
Rgam2 = Sqrt(deltar2((1))^2 + deltar2((2))^2);

Ec2 = (e/(4 Pi*epsilon)) (deltar2/Rgam2 - betas)/(gamma^2 Rgam2^2 (1 - (deltar2/Rgam2).betas)^3);
Erad2 = (e/(4 Pi*epsilon)) Cross(deltar2/Rgam2,Cross(deltar2/Rgam2 - betas, betaDot))/(c*Rgam2*(1 - (deltar2/Rgam2).betas)^3);

sumElong = (Ec2((1)) + Erad2((1)));
sumEtran = (Ec2((2)) + Erad2((2)));

(*Numerical Functions*)

ZPRIME(NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, s_?NumericQ, sprime_?NumericQ) := zprime /.FindRoot(s - sprime == (Sqrt(gamma^2 + KK^2) (EllipticE(kw*(NN + 10000)*lw,KK^2/(gamma^2 + KK^2)) - EllipticE(kw zprime, KK^2/(gamma^2 + KK^2))))/(gamma kw) -beta Sqrt(((NN + 10000)*lw - zprime)^2 + (x - xprime + (KK Sin(kw *(NN + 10000)*lw))/(gamma kw) - (KK Sin(kw zprime))/(gamma kw))^2), {zprime, 0})

Kernel = coeff re/gamma (sumElong*betao((1)) + sumEtran*betao((2)))/.{zprime -> ZPRIME(NN, x, xprime, s, sprime)};

FNxprimesprime(NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, s_?NumericQ, sprime_?NumericQ):= Kernel

FNsprime(NN_?NumericQ, x_?NumericQ, s_?NumericQ, sprime_?NumericQ) :=NIntegrate(FNxprimesprime(NN, x, xprime, s, sprime), {xprime, -300/10^6, 300/10^6})

FN(NN_?NumericQ,x_?NumericQ, s_?NumericQ) := NIntegrate(FNsprime(NN,x, s, sprime), {sprime,-10^-4, s})

lst1 = Table({ss, FN(0,0, ss), PrecisionGoal -> 5) // Quiet}, {ss, -10^-5, 10^-5, 10^-6})