Tag Info

New answers tagged

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}]]; ...


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[ ...


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 ...


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, ...


0

Modifying your second definition of FT to FT[kx_, ky_, z_] := 1/(2*Pi)* NIntegrate[ f[x, y, z]*Exp[-I*kx*x - I*ky*y], {x, -∞, ∞}, {y, -∞, ∞}, MinRecursion -> 4] // Chop solves your problem.


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 ...


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]; ...


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}] // ...


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: ...


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 ...


13

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] := ...


3

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" ...


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 ...


0

You can use the finite element method with the method of lines as @toadatrix suggested, but for the FEM method to work, you need to do a little more. The Neumann boundary conditions need to be specified using NeumannValue. h[x_] := x*(30 - x)/900; op = D[u[t, x], t] - D[u[t, x], x, x]; begin = 0; end = 30; bc = {u[0, x] == 100*h[x]}; neumann = ...


0

I think this might have been a bug. At least as of V10.0.0, the following works without incident: simplesys = {r'[t] == Piecewise[{{1, 0 <= t <= 10}, {0, 10 <= t <= 20}}, 0], r[0] == 0}; {state} = NDSolve`ProcessEquations[simplesys, r, t]; Here we reinitialize three times and plot the results of each: Module[{newstate, sol}, ...


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}]


8

I think there's a bug in the internal function NDSolve`SPRKDump`CheckSeparability 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 ...


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] ...


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}]


0

In the absence of specific definitions for the matrices in the question, I simply made up my own. Here is a way you can make it work: J = {{1, 2}, {2, 1}}; T = 1; γi = {{1, 0}, {0, 0}}; g = Identity; Clear[y, Y] y = Function[t, ##] &[Array[Y[##][t] &, Dimensions[J]]]; eqn1 = Thread /@ Thread[y'[t] + (J\[Transpose] + ...


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 ...


5

With two minor changes to the code bm = ToBoundaryMesh["Coordinates" -> coord, "BoundaryElements" -> {LineElement[incidents]}]; nds = NDSolveValue[{Inactive[Div][sigma*Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0, DirichletCondition[u[x, y] == 0, y == 0], DirichletCondition[u[x, y] == 1, y == 1]}, u, {x, y} ∈ mesh]; the ...


1

[Background: I added this answer because this question was cited in What are the most common pitfalls awaiting new users? as an example of when NumericQ was required. I took the integrand to be a simple example of a more general problem with NIntegrate, but on reading the comments, it seems to be narrowly the exact integrand of interest to the OP. I ...


3

The error messages come from Integrate/NIntegrate calls on functions with symbolic (non-numerical) parameters a, b, c. These calls come before FindRoot substitutes numerical values for these parameters. (Taking N of Integrate here is effectively the same as calling NIntegrate.) In this case the best approach it seems to me is to do the integral first, ...


5

The the function y is the integral of a logarithmic singularity, so it is relatively easy to understand. The function JacobiCN[s, 7/10] is relatively flat near s == EllipticK[7/10], so that the integral of its reciprocal can be approximated by a logarithm: y[s]/JacobiCN[s, 7/10] ~= y0 / (jp (s - EllipticK[7/10])) where y0 ~= y[s] at s == EllipticK[7/10] ...


5

My question is: how to set that NDSolve will not save whole InterpolationFunction for the result? There is actually a very simple way to do this: instead of specifying a list of functions in the second argument, specify an empty list instead. This now begs the question of how one can obtain results. The solution is to use the event location ...


12

WhenEvent is working. Try WhenEvent[x[t] < 0, Print[t]; x[t] -> 0] to see that every crossing is detected. The problem is that it only detects crossings. So changing x[t] -> 0 does not reset the event. At the next step x[t] becomes negative and no event is detected. (This is how it is supposed to behave.) The way to deal with this is to use ...



Top 50 recent answers are included