Tag Info

46

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

41

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["NDSolveFEM"] mesh = ToElementMesh[FullRegion[2], {{0, 5}, {0, 1}}, "MeshOrder" -> 1, ...

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

24

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

23

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

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

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

22

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

18

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

18

You can use the EventLocator method of NDSolve. Needs["DifferentialEquationsInterpolatingFunctionAnatomy"]; 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

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

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

17

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

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

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

17

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

16

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

16

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

16

NIntegrate performs a certain symbolic processing of the integrand to detect discontinuities, singularities, to determine the method to choose and so on. If you know the integrand pretty well, the way to reduce the overhead is to set the method explicitly, set its SymbolicProcessing suboption to 0 (to allow to time spent on the preprocessing), and to add ...

16

As it turns out, the designers of NDSolve[] have precisely anticipated this sort of use; this is where you can use the NDSolveStateData framework. To use acl's example: (* prepare PDE *) state = First[NDSolve`ProcessEquations[{D[u[t, x], t] == D[u[t, x], x, x], u[0, x] == 0, u[t, 0] == Sin[t], u[t, 5] == 0}, u, t, {x, 0, 5}]]; (* go up to t = 2 *) ...

16

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

16

When n is large it's much faster to operate on a 3 x n array than to process each of the n 3-vectors separately. This is one of the standard "tricks" to speed things up. n = 10^6; (* Isn't that easier to read than 1000000 ? *) AbsoluteTiming @ N[ Total@UnitStep[ 1. - Norm/@RandomReal[1,{n,3}] ]/n ] (* {4.555842, 0.524302} *) AbsoluteTiming @ N[ ...

15

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

15

You can get the curve in polynomial implicit form as below. poly = GroebnerBasis[{x^2 - ct, y^2 - st, ct^2 + st^2 - 1}, {x, y}, {ct, st}][[1]] (* Out[290]= -1 + x^4 + y^4 *) To get the area, integrate the characteristic function for the interior of the region. That that's where the polynomial is nonpositive (just notice that it is negative at the ...

15

According to the Mathematica documentation on this page: Here is how to define a 5(4) pair of Dormand and Prince coefficients [DP80]. This is currently the method used by ode45 in MATLAB. DOPRIamat = { {1/5}, {3/40, 9/40}, {44/45, -56/15, 32/9}, {19372/6561, -25360/2187, 64448/6561, -212/729}, {9017/3168, -355/33, 46732/5247, 49/176, ...

15

Here's my attempt. To get the matrix representing the Laplacian I use LaplacianFilter on an array of symbols and CoefficientArrays to extract the coefficients. n = 200; shape = ArrayPad[ConstantArray[0, {n/2, n/2}], {{0, n/2}, {0, n/2}}, 1]; shapeVector = Flatten @ Position[Flatten @ shape, 1]; symbolArray = Array[x, {n, n}]; symbolLaplacian = ...

14

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

14

NDSolve has a slew of options that allow you to control the method. You can find the standard reference here. There, we learn how to access Euler's method using NDSolve: Clear[x]; x = x /. First[ NDSolve[{x'[t] == 0.5*x[t] - 0.04*(x[t])^2, x[0] == 1}, x, {t, 0, 10}, StartingStepSize -> 1, Method -> {"FixedStep", Method -> "ExplicitEuler"}] ]; ...

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