# Tag Info

## Hot answers tagged numerical-integration

14

The memory leak in NIntegrate is a bug and has been fixed as of version 10.2.0. Earlier versions would lose around 720 bytes per evaluation for this example, which could not be recovered without restarting the kernel. ClearSystemCache[] should be used to make sure the memory is released. Using version 10.2: NI[z_?NumericQ, b0_?NumericQ] := ...

8

I think there's a bug in the internal function NDSolveSPRKDumpCheckSeparability that leads NDSolve to conclude that the system is not separable. I think you should report it and see if WRI can verify it (they would probably appreciate a link to this Q&A). It's a fair amount of work to track it down, and there is a lot of nearly unreadable stuff to ...

6

This problem can be simplified enormously by performing the integral over ϕ symbolically. To do so, first consider the ϕ-dependent exponential in the integrand. Exp[I k (x Sin[θ] Cos[ϕ] + y Sin[θ] Sin[ϕ] + z Cos[θ])] It can be transformed and simplified as follows. expt = Exp[I k ( x Sin[θ] Cos[ϕ] + y Sin[θ] Sin[ϕ] + z Cos[θ])] /. {x -> r ...

6

Using polar coordinates r and f, the region of integration is given by { 0<= r <=2/Cos[f], 0<= f <= 2 \[Pi] } We can then proceed as follows. First integration with PrincipalValue g = 8 Integrate[r/(1 - r^2), {r, 0, 2/Cos[f]}, Assumptions -> 0 < f < \[Pi]/4, PrincipalValue -> True] (* Out[451]= -4 Log[-1 + 4 Sec[f]^2] *) ...

5

Your problem is partly that NIntegrate calls a function that at the given point in time is not full numerical (the first expression in fn). A way around that is to define the functions in such a way that hey will only evaluate for purely numerical values: Clear[m1, m2, fn, gn]; m1[x_?NumericQ] := (2 (x - 1))/x; m2[x_?NumericQ, y_?NumericQ] := (x (2 - x ...

5

Compile f and use a memo-ized version of it Since it seems like NIntegrate decides to symbolically evaluate its argument first, I thought I'd force it not to by compiling the function f. This seems to make a significant difference: Clear[f, f1, g] g[x_] = Nest[f[x] + 1./# &, f[x], 500]; f1 = Compile[{x}, Sum[1/100 Erfc[-(x^2/k)], {k, 100}]]; ...

4

Integrate can be used to handle the DiracDelta, Integrate[Integrand[p, Q2, ν, θ], {p, 0, Infinity}, Assumptions -> p3zero2[Q2, ν, θ] ∈ Reals && p3zero1[Q2, ν, θ] ∈ Reals] (* ((HeavisideTheta[p3zero1[Q2, ν, θ]]*p3zero1[Q2, ν, θ]^2)/ ((M - 2*E3[p3zero1[Q2, ν, θ]])^2*E3[p3zero1[Q2, ν, θ]]* E4[p3zero1[Q2, ν, θ], Q2, ν, θ]*Derivative[1, 0, ...

4

To be honest, the code you found isn't a good example of coding in Mathematica. I think the following 3 lines are enough for you: ListStreamPlot[Transpose[{u, v}, {3, 1, 2}], DataRange -> {{0, 1}, {0, 1}}] curl = Most /@ Differences@v - Most@Differences[u, {0, 1}]; ListContourPlot@LowpassFilter[curl, 1] You may also want to try ListDensityPlot: ...

4

Update 3 As ilian correctly states in the comments, without NIntegrate we get the same picture as shown in the "Update 2" section. It proves that memory grow reflects temporary memory allocation by Table, and ClearSystemCache[] does completely eliminate the memory leak. We can also switch off the caching completely: SetSystemOptions[ "CacheOptions" ...

4

While the unknown g might or might not be amenable to this approach, it's really fast on g == 1: Expectation[ 1, {x1, x2, x3, x4} \[Distributed] MultinormalDistribution[muvec, sigmat]] // RepeatedTiming (23000 times faster than NIntegrate.) Then there's also NExpectation to try, too. Ever since I came across Guess who it is's use of the ...

3

Because the integrand is highly localized at the origin; e.g., Log[10, npdf[{6, 0, 0, 0}]] // N (* -9.16177 *) limiting the range of integration reduces run time by a factor of three and moderately improves accuracy. muvec = ConstantArray[1/10, 4]; sigmat = IdentityMatrix[4]; npdf[x_] := PDF[MultinormalDistribution[muvec, sigmat], x]; ...

3

I do not know what range of $s$ you had in mind, so I arbitrarily chose one below. You should be able to change it at will. This is an example in hopes that you will study this further and get to the point where you can come up with such code yourself. Please take a look at the excellent links provided by Arnoud as well to get you started. Clear[phi] ...

3

Setting WorkingPrecision -> 5 in the gn integral gives you a reasonable convergence time. At the expense of some more computation time you can check that methods DifferentialEvolution and SimulatedAnnealing both return the sme result given here up to four decimal places. m1[x_] := 2 (x - 1)/x; m2[x_, y_] := x (2 - x y)/(2 (x - 1) y); fn[x_?NumericQ, ...

2

You can integrate as follows. Integrate[BesselI[-nu, k*x]/x, {x, r1, r}, Assumptions :> {k \[Element] Complexes, r1 \[Element] Reals, r \[Element] Reals, nu \[Element] Reals, r1 > 0, r > r1} The result is a complicated expression in terms of Gamma and HypergeometricPFQRegularized functions. Nevertheless, it can be ...

2

It is quite tricky! Piecewise[ ] functions work only with "EquationSimplification" -> "Residual"... I'll try to dig up why simplesys = {r'[t] == Piecewise[{{1, 0 <= t <= 10}, {0, 10 <= t <= 20}}, 0], r[0] == 0}; state = First@ NDSolve`ProcessEquations[simplesys, r, t, Method -> {"EquationSimplification" ...

2

Here you have a way without using neither UnitStep nor PieceWise that improves the performance by 75% wrt your code. It computes 200 functions in a very reasonable time for your toy example. The main trick is to use a numeric (black box) function to be able to take Part[... ] inside it. n = 200; af = Array[f, n]; taf[t_] := Through[af[t]] bb[i_?IntegerQ, ...

2

ClearAll[δ, f]; f := 2 Δ^2 Sin[(A + B)/2]^2 Re@NIntegrate[ E^(-I t (2 ϵ + F + G)/Λ) E^(-Z/Λ (t ArcTan@t- Log[1 + t^2]/2)), {t, 0, ∞}] Plot[f, {δ, 0, 100}]

2

The reason for this error is that NIntegrate uses fixed precision when computing the integration ranges, while EllipticK needs to raise the precision internally to obtain a good result. N[EllipticK[7/10], 20] (* 2.0753631352924691439 *) Block[{$MinPrecision =$MaxPrecision = 20}, N[EllipticK[7/10], 20]] (* Divide::infy: Infinite expression ...

2

Clear[integrand] With your definition that uses an unnecessary list and then takes a part integrand[t_?NumericQ] := {1000/Sqrt[t^2 - tMin0^2], 0}; integrand[t][[1]] t This occurs because integrand cannot evaluate with a symbolic input so the first part of the unevaluated expression is the function's argument. Clear[integrand] tMin0 = 1/10; tMax0 ...

2

The following part is incorrect, and redundant anyway: \[Psi]= \[Integral]u dx= \[Integral]v dy \[Omega]= \[PartialD]v/dx-\[PartialD]u/dy Otherwise, the code you copied has only one major problem, in the definition of n. If you try to execute it, you get Do::iterb: Iterator {j,2,n-1} does not have appropriate bounds. >> which suggests where the ...

1

The first thing I did was to rationalizing all calculations, starting with the defintion of σ and minroot. This stops the Solve::ratnz messages. I also made some other improvements to minroot. σ = 6/10; minroot[gg_?NumericQ, bb_?NumericQ] := Module[{b, g, rts, r}, b = Rationalize[bb, 0]; g = Rationalize[gg, 0]; rts = r /. Solve[ ...

1

At the moment the LUM function as written can't be used because it depends upon a number of global variables (h, me, k, Tcmd ...) that you have not supplied. However, in order to produce a list of {x2, LUM[x1,x2]} pairs this can be done as follows. Scalar Assume x2List is the list of values that you want to evaluate and you want to evaluate LUM at a ...

1

About 30X faster integration, no errors... Clear[x, x1, x2, x3, x4] muvec = ConstantArray[1/10, 4]; sigmat = IdentityMatrix[4]; pdf = PDF[MultinormalDistribution[muvec, sigmat], {x1, x2, x3, x4}] // FullSimplify Integrate[ 1*pdf, {x1, -Infinity, Infinity}, {x2, -Infinity, Infinity}, {x3, -Infinity, Infinity}, {x4, -Infinity, Infinity}] // ...

1

This is your function: phi[s_, H_, B_] := NIntegrate[(1 + B Sin[2 H t])/Sqrt[ 1 + B^2 + 2 B Sin[2 H t]], {t, 0, s}] 1/(2 H) Sqrt[ 1 + B^2 + 2 B Sin[2 H s]] But if you want to plot it you have to give the value for s. Plot[{phi[s, 0.5, 1], phi[s, 1, 1], phi[s, 1.5, 1], phi[s, 2, 1]}, {s, 0, 16}]

Only top voted, non community-wiki answers of a minimum length are eligible