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


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


15

It turns out ListSurfacePlot3D does a terribly poor job of approximating the surface in the OP, otherwise one will just apply DiscretizeGraphics to the output obtained from ListSurfacePlot3D and be done with it. But since that's not applicable here, we present an approach that uses alpha shapes to approximate the shape of the given point set by tuning a ...


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


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

NDSolve currently can't handle this kind of differential equation, LaplaceTransform is your friend. Since in this case inverse Laplace transform can't be done analytically by InverseLaplaceTransform, you need the help of numerical Laplace inversion package in addition: eq = {0.01 - 6.25 x[t] + (1.2 Integrate[x'[t - τ]/Sqrt[τ], {τ, 0, t}])/10^7 == 16 ...


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


12

$Version (* Out[228]= "8.0 for Microsoft Windows (64-bit) (October 7, 2011)" *) The result provided by Mathematica is correct: Integrate[a/(a^2 + Sin[t]^2), {t, 0, 2 π}] (* Out[213]= (2 π)/(Sqrt[1 + 1/a^2] a) *) Now the same procedure as "always" which "explains" the zero result. The indefinite integral is Integrate[a/(a^2 + Sin[t]^2), t] (* ...


12

For what it's worth, you can monitor the progress of Integrate and NIntegrate. I'm not sure how helpful the tools below are, but I feel it is probably worth mentioning them. Integrate Internal`Integrate`debugSwitch If you set Internal`Integrate`debugSwitch to the magic number 10, it will print its (major) steps in searching for the answer. For instance: ...


12

If you're running Mathematica 10 or above, there's a dead simple method using the new Area function: newimg = img /. Ball[x___, 0.6] -> Sphere[x, 0.5] Total[Map[Area[RegionIntersection[#, Cuboid[{-1, -1, -1}, {1, 1, 1}]]] &, newimg]] (* 25.1327 *) Note a few things here: Mathematica considers Ball to be a 3D object; if we want its surface area ...


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

Working solution One can manually implement the shooting method with ParametricNDSolveValue and FindRoot: psol = ParametricNDSolveValue[{k f[z] + f'[z] == 0, f[-1] == Exp[2], WhenEvent[z == 0, f[z] -> r0 f[z]]}, f, {z, -1, 1}, {k, r0}]; k0 = k /. FindRoot[psol[k, 0.5][1] == 1, {k, 1}] (* 0.653426 *) Plot[psol[k0, 0.5][z], {z, -1, 1}] ...


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


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


8

One can use Picard-type iteration to get the solution: Using an approximation to x'[t] (in the integral), we can integrate the ODE to obtain a new approximation. Remarkably, it converges in just two steps. My original thought was to step through the integration using the tools from tutorial/NDSolveStateData to build an interpolation of x'[t] at each step ...


8

It looks like classic catastrophic round-off error. (Look at those exponents on $e$!). {sol} = DSolve[{x'[t] == y[t]/100, y'[t] == -100 x[t] - 100 y[t] + 2020, x[0] == 0, y[0] == 20}, {x, y}, t]; Now consider y[9.] vs. y[9] and y[9.`20]: y[9.] /. sol N::meprec: Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating ...


8

You can make a change of variable to solve the problem. Here I'll use dChange for this task: r0 = 0.5; eqn = {k[z] f[z] + f'[z] == 0, k'[z] == 0, f[-1] == Exp[2], f[1] == 1}; c = Piecewise[{{r0, z > 0}}, 1]; neweqn = dChange[eqn, f[z] == c g[z]]; {solg, solk} = NDSolveValue[neweqn, {g, k}, {z, -1, 1}]; Plot[c solg[z], {z, -1, 1}]


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*)



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