Hot answers tagged

35

General comments First, if you plan to use multi-dimensional integrals it is better to test with multi-dimensional integrals not with one dimensional ones. One might think that the test in the question is an appropriate one if multi-dimensional integration is done by the integrator in a recursive manner. This seems to be case for scipy.integrate.nquad (see ...


24

The simplest way to make new NIntegrate algorithms is by user defined integration rules. Below are given examples using a simple rule (the Simpson rule) and how NIntegrate's framework can utilize the new rule implementations with its algorithms. (Adaptive, symbolic processing, and singularity handling algorithms are seamlessly applied.) Basic 1D rule ...


20

This question comes up often enough. See this discussion at community.wolfram.com : Integration method used in NIntegrate , and the notebook Finding the applied NIntegrate methods attached to my second response in the discussion. That notebook contains examples of usage of the undocumented function NIntegrateSamplingPoints and NIntegrate's option ...


19

I'll preface this answer first with a complaint: NExpectation[] and NProbability[] are not sufficiently resilient obviously adjustable. Ideally, these two functions are an "interface" to NIntegrate[], allowing the user to formulate his expression purely in distributional terms. Unfortunately, when one hits cases like this, the things one might usually ...


19

Motivation (for a new semi-symbolic integration strategy) Consider the following integral, which cannot be done neigther by Integrate: Integrate[BesselJ[y, x^3], {x, 0, ∞}, {y, 0, 1}] (* Integrate[If[Re[y] > -(1/3), Gamma[1/6 + y/2]/(3*2^(2/3)*Gamma[5/6 + y/2]), Integrate[BesselJ[y, x^3], {x, 0, Infinity}, Assumptions -> Re[y] <= -(1/3)]...


15

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] := NIntegrate[E^-...


15

Introduction I think there are several questions on this site about ODEs of the form $$(x-a)^2 u''(x) = F(x,u,u')$$ with an initial condition at $x=a$. There is no general guarantee that solutions exist over an interval $(a,b]$, but sometimes it is possible as in this case. Outline We transform the equation $u''(x) = F(x,u,u')$ over the infinite interval ...


14

In my experience FindRoot works best for such problems: In[1]:= fun[a_?NumericQ] := NIntegrate[(x/(Exp[x] - 1)), {x, 0, a}] In[2]:= FindRoot[1 - (2/a)*((1/a)*fun[a] + (a/6) - 1) == 0.7, {a, 0.1}] Out[2]= {a -> 58.3073}


13

As a minor addition to J.M.'s excellent answer, If something breaks while evaluating NExpectation[] or NProbability[], a number of the things that otherwise can be adjusted in NIntegrate[] aren't there. Options can be passed to NIntegrate, for example try something like Table[NExpectation[X, X \[Distributed] JohnsonDistribution["SB", γ, δ, 0, 1], ...


13

The main problem is that your pos is not seen as a 3D vector. The cross product is therefore interpreted as a scalar: q*Cross[D[pos[t], t], b] when adding this to the vector q.e this 'scalar' term is added to each of the vector components: q*e + q*Cross[D[pos[t], t], b] This won't work, instead do: b = {1, 0, 0}; e = {0, 0, 1}; q = 1; m = 1; ...


13

Some explanations first The substitution in the question introduces the reduced wave function $u(r)$ by solving the original radial equation in polar coordinates, $$-\frac{1}{2}\left(R''(r)+\frac {1}{r}R'(r)\right) - \frac{1}{r}R(r) + \frac {m^2}{2r^2}R(r) = E R(r)$$ using the ansatz $$R(r)\equiv \frac{1}{\sqrt{r}}u(r)$$ The apparently divergent ...


12

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


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


11

Because of the comment by @ciao I think it is a good idea to give a solution/answer using MonteCarloRule and PointGenerator. PointGenerator is an object that expects a sub-value function definition for "MakePoint" with the signature: "MakePoint"[dim_, totalNumberOfPoints_, i_, wprec_] where the arguments are: dim - integration dimension, ...


11

s = ParametricNDSolveValue[{x'[t] == -y[t] + x[t]*Log[x[t]], y'[t] == x[t] + y[t]*Log[x[t]], x[0] == x0, y[0] == 0}, {x, y}, {t, 1}, x0] f[x0_, t_] := Through[Through[s@x0]@t] pts = Table[f[x0, t], {x0, 1, 2, .2}, {t, 0, 1, .1}]; Show[Graphics[{Green, Arrow /@ pts, Black, Point /@ pts}, ...


10

(Update: I forgot to copy some of the code) For values of z = x beyond the critical value (just below z == 7500), there are three real equilibria, two stable spirals and one (unstable) saddle. eq[x_] := Module[{q = x}, w0 = 7000; G0 = 50; Q = 14000; hU = 0.6*G0; w00 = w0 - q + Q*Abs[A0[t]]^2; At = 10^-6; G1 = 1.0 G0; tmax = 4; e1 ...


10

EDIT #2 My error was useful. It brought me to the conclusion that the difficulties in solving the PDE of the OP are due to the drift term $$\frac{\partial (x u(x,t))}{\partial x}$$ If the drift term is included, many boundary problems are ill defined. It turns out that there are cases where mathematically there is only a trivial solution u = 0 but ...


10

The following is not particularly fast and could be more accurate but does make progress toward the goals set in the Question. To begin, consider the analytical properties of f, as defined in the Question. By observation, it has branch points at {I Sqrt[6/5] ξ, I Sqrt[2/5] ξ} and their conjugates. Poles are obtained by p /. Simplify[Solve[...


9

Mathematica seems to split the integrand component, E^(-((-m + Log[x])^2/(2 s^2))) into E^(-((m^2 + Log[x]^2)/(2 s^2))) times the sort-of "coefficient" E^((m Log[x])/s^2) (* == x^(m/s^2) *) in order to calculate the integral in terms of Meijer $G$. For reasons that are obscure to me, it seems to want the coefficient of m in the exponent to be ...


9

First, in general, I would advise you not to trust numerical algorithms. If there are doubts about the outcomes then solve the same problem with different (numerical or not) methods and see do their results agree. For the integral in the question I assume you can evaluate it with several different invocations of the Monte Carlo method and compare the ...


9

This is the best I can come up with, I'm very interested to see if anyone else has a better solution. The idea here is to just run through values of $t$, and do a DFT on $$E(t+\frac{\tau}{2}) E ^*(t-\frac{\tau}{2})$$ So I set up the time/frequency resolution for my DFT, using a dt value I know gives a broad enough spectrum, dt = 0.025; num = 2^14; df = ...


9

Using the undocumented IntegrationMonitor: {val, {vals}} = Reap@NIntegrate[x y^2, {x, 0, 1}, {y, 0, 1}, PrecisionGoal -> 2, Method -> "MonteCarlo", IntegrationMonitor :> ((Sow[Total@Through[#["Integral"]]]) &)] (* {0.165268, {{0.186623, 0.172189, 0.168129, 0.166339, 0.165429, 0.16988, 0.173145, 0.171675, 0.173355, 0.177199, 0....


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


8

This should work: Clear[f]; f[a_?NumericQ, b_?NumericQ] := NIntegrate[Sqrt[(Cos[t] - a)^2 + b^2], {t, 0, Pi}] and then add //N at the end of the definition of g[a,b] g[a_, b_] := Derivative[1, 0][f][a, b]//N g[1,1] (*1.80525*)


8

I'm not confident of this answer, but it seems worth presenting: NIntegrate[2^3/((x - s)^2 + (y - t)^2 + (z - u)^2), {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, {s, 0, x}, {t, 0, y}, {u, 0, z}, AccuracyGoal -> 16] // Quiet (* 5.63378 *) The basic insight is that each "paired" integrations (e.g., $x$ and $s$; $y$ and $t$; $z$ and $u$) is over the unit square. ...


8

The trick is to use NDSolve instead of NIntegrate and thus in effect obtain a numerical antiderivative that can be evaluated fast at different points. NIntegrate will only do definite integration, so it needs to be run each time the integration bounds are changed. This is very slow, as you noticed. NDSolve will only need to be run once. Slow way (what ...


8

Having put in some time trying to see what's going, I've found a few things, but I don't have a perfectly clear picture. I believe the issue is with the large InterpolatingFunction in the integrand and not with NIntegrate per se. The time it takes for NIntegrate to set up the integration is much longer in V10, but the integration itself runs in about the ...


8

Since f[t] yields in effect {x[t], y[t]} and the area under a curve (with y[t] > 0) is given by the integral of y[t] x'[t] (i.e., $\int y\;dx$), then the following should work, assuming x[t] is increasing. dA[f_, t_?NumericQ] := Last[f[t]] First[f'[t]]; NIntegrate[dA[f, t], {t, 0, 1}] (* 150000. *) One can extend this to other forms, e.g. to closed ...


8

Alternative method: b = {1, 0, 0}; e = {0, 0, 1}; q = 1; m = 1; sol = NDSolve[{e + Cross[pos'[t], b] == m/q pos''[t], pos[0] == {0, 0, 0}, pos'[0] == {0, 0, 0}}, pos, {t, 0, 10}, Method -> {"EquationSimplification" -> "Residual"}]; ParametricPlot3D[pos[t] /. sol, {t, 0, 10}, PlotRange -> All]


8

We calculate the symbolic result of the integral for the two values of the parameter δ. Considering, as I did before in a related problem (Strange result of parameter-dependent definite integral), the exponential form instead of the cosine I obtain after expanding the exponential function containig δ and integrating term by term the following result (where ...



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