Hot answers tagged

53

Control the Precision and Accuracy of Numerical Results This is an excellent question. Of course everyone could claim highest accuracy for her product. To deal with this situation there exist benchmarks to test for accuracy. One such benchmark is from NIST. This specific benchmark deals with the accuracy of statistical software for instance. The NIST ...


47

The only reason I am attempting to answer this is to perhaps get a Reversal badge. There you go... We will go slowly and this answer is the basis for what comes next. Let's start with two dimensions. You'll see why. We create a rectangular region: Needs["NDSolve`FEM`"] mesh = ToElementMesh[FullRegion[2], {{0, 5}, {0, 1}}, "MeshOrder" -> 1, "...


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


30

Solving 1D and 2D complex Schroedinger wave equations with NDSolve I do not agree with you when you write: I know the NDSolve is not magic... My opinion is that NDSolve is one of the most complex functionality I've met so far in the Mathematica environment, with its millions of options and special function this is a real complex thing and it is hard ...


29

Let me show how to roll your own numerical solution to a non-linear integral equation using a collocation method. It's fun! This will involve two approximations. First, we will approximate the function B[x] by its values at n particular points in the range {x, 0, 1}. The integral over x will be replaced by a weighted sum over n, i.e., a quadrature rule. ...


27

Time-dependent case in the time-dependent case, $[H(t),H(t')]\neq0$ in general and we need to time-order, ie, the operator taking a state from $t=0$ to $t=\tau$ is $U(0,\tau)=\mathcal{T}\exp(-i\int_0^\tau dt\, H(t))$ with $\mathcal{T}$ the time-ordering operator. In practice we just split the time interval into lots of small pieces (basically using the ...


25

There is an (undocumented?) feature of NDSolve which is handy for exactly this purpose: You can add more than just the start and end of the integration interval and enforce that these points will be met. The result is like you would run NDSolve on each of the corresponding intervals with the starting conditions given by the end point of the previous interval....


25

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


23

Edit of July 10, 2014 As of V10, this equation can now be solved with a single, simple call to NDSolve: y = NDSolveValue[{ r D[y[r, z], z, z] + D[y[r, z], r] + r D[y[r, z], r, r] == r y[r, z], y[1, z] == 1, y[r, 1] == 1 }, y, {r, 0, 1}, {z, 0, 1}]; ContourPlot[y[r, z], {r, 0, 1}, {z, 0, 1}, ColorFunction -> "TemperatureMap", ...


23

The programming style you are using is not very fitting for Mathematica. Here's a better way (shorter, much faster): n = 1000000; (* number of points to use *) octantVolume = N[ Total@UnitStep[1 - Norm /@ RandomReal[1, {n, 3}]]/n ] The reason why you get the error you mention is that for some x, y, the expression 1 - x^2 - y^2 is negative, thus its ...


22

You can always separate your inner integrals, convert them to functions and use in NIntegrate: i1[z_?NumericQ] := i1[z] = NIntegrate[-y, {y, 0, z}] i2[x_?NumericQ] := i2[x] = NIntegrate[Exp[i1[z]], {z, -∞, x}] NIntegrate[x i2[x], {x, -5., 5}] (* 30.0795 *)


21

I think it's worth pointing out that the problem can be solved "straightforwardly" (i.e., really using only NDSolve) once you know the options that Stefan used in ProcessEquations (which I upvoted because those options are the main ingredient): Below I show the original problem of a Gaussian wave packet with no initial momentum, and then a modified case ...


20

You can express your integral in terms of a differential equation and use NDSolve. Since NDSolve builds up the solution as it goes, this is typically much faster. Clear[y]; y[x_] = y[x] /. First[ NDSolve[{y'[x] == Sin[x], y[0] == 0}, y[x], {x, 0, 10}] ]; t = AbsoluteTime[]; Plot[y[t], {t, 0, 10}] AbsoluteTime[] - t


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


20

After a lengthy study (I'm using version 8) I conclude that there is a bug in Mathematica in the Integrate function when applied to a Sqrt integrand. Ok. let's go (some patience is required because of the long text) Let us define the functions corresponding to your integrals. Remark: because of the relation $1 + cos(2x) = 2 cos^2(x)$ the two forms of ...


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


18

There is the function NFourierTransform[] (as well as NInverseFourierTransform[]) implemented in the package FourierSeries`. The function, as with the related kernel functions, takes a FourierParameters option so you can adjust computations to your preferred normalization as needed. For your specific normalization, you apparently want the setting ...


18

My variant of Szabolcs code. It doesn't need an extra package: sol = First[ NDSolve[eqns, {a, b}, {t, 0, 1000}, Method -> {"EventLocator", "Event" -> Abs[a'[t]] +Abs[b'[t]] < 10^-5, "EventAction" :> Throw[end = t, "StopIntegration"]}]]; Plot[Evaluate[{a[t], b[t]} /. sol], {t, 0, end}] As you can see it makes use of the "...


18

You can use the EventLocator method of NDSolve. Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]; eqns = {Derivative[1][a][t] == -a[t] - 0.2` a[t]^2 + 2.1` b[t], Derivative[1][b][t] == a[t] + 0.1` a[t]^2 - 1.1` b[t], a[0] == 0.5`, b[0] == 0.5`}; sol = First@ NDSolve[eqns, {a, b}, {t, 0, 1000}, Method -> {"EventLocator", "...


18

Some frames from my version of the animation: Here's the code I used: orbit[posStart_?VectorQ, derStart_?VectorQ] := Block[{c = -Rationalize[6.672*^-11*7*^17], x, y, z, t}, {x, y, z} /. First @ NDSolve[ Join[Thread[{x''[t], y''[t], z''[t]} == c {x[t], y[t], z[t]}/Norm[{x[t], y[t], z[t]}]^3], ...


18

Method of random number generation is also significant: Default: n = 10^6; AbsoluteTiming[N@Mean@UnitStep[1. - Total[RandomReal[1, {3, n}]^2]] - π/6] {0.197896, 0.000649224} Niederreiter low-discrepancy sequence (see "methods" here): SeedRandom[Method -> {"MKL", Method -> {"Niederreiter", "Dimension" -> 3}}]; AbsoluteTiming[N@Mean@UnitStep[...


17

Looks like a bug that ran off in the development version of Mathematica. In[1]:= Integrate[Exp[I Cos[b - c]] Cos[b], {b, 0, 2 Pi}]//InputForm Out[1]//InputForm= (2*I)*Pi*BesselJ[1, 1]*Cos[c] In[2]:= %/.c->.5 Out[2]= 0. + 2.42645 I


17

Generally speaking, you can recognize a list because it'll have List as its Head. For example: Head[{1,2,3}] will return List. For your example conditional where you want to change what you do based on the Head of the resulting expression, you can use Switch, such as in: Switch[result, _List, what you want to do with a list, _, what you ...


17

Mathematica is an incredible tool for checking conjectures and making sketches. I'm going to demonstrate it below. Let's start with checking that in case when $R_0$ (I replaced it with $R$) is a polynomial, the solution of this Volterra equation reduces to linear ODE. Lets take some derivatives of the equation: ClearAll[P, R, s, t]; eqn = P[t] == R[t] + ...


17

You can modify the global system variable $Assumptions, to get the effect you want: $Assumptions = aa[t] > 0 Then Integrate[D[yy[x, t], t]^2, {x, 0, 18}] 10.1601 Derivative[1][aa][t]^2 This may, however, be somewhat error-prone. Here is how I'd do this with local environments. This is a generator for a local environment: createEnvironment[...


17

There's absolutely no need to load a package if all you want to do is simple Gauss-Legendre quadrature: GaussLegendreQuadrature[f_, {x_, a_, b_}, n_Integer: 10, prec_: MachinePrecision] := Module[{nodes, weights}, {nodes, weights} = Most[NIntegrate`GaussRuleData[n, prec]]; (b - a) weights.Map[Function[x, f], Rescale[nodes, {0, 1}, {a, b}]...


17

If you insert a bunch of commands like "Print@First@AbsoluteTiming[...]" in the middle of your function you will see that literally all time is spent on the last line with the Sum[NIntegrate[...]]. To solve your problem insert the following commands into NIntegrate: Method -> {Automatic, "SymbolicProcessing" -> 0} Caution: I have personally ...



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