# Tag Info

42

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

26

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

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

22

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

20

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

19

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

17

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

17

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

17

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

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

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

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

15

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

15

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

15

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

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

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

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

13

I think you intended to use {li, 200, 800} instead of {li, 800, 200}. If you do so, then you could visualize the result : ListLinePlot@dnFpoints Moreover I would rather define daF in the following form : daF[l_]:= 500 * 0.28 Exp[-((l - 500)/90)^2] c = 3 10^8; Edit Instead of using Table of dnFpoints I add an alternative method for calculation of ...

13

This question is somewhat subjective, but here's my take on it: The reason the precise methods are mentioned in papers is to make results reproducible. One has to draw a line when it comes to describing methods. Will you mention what method you used to add or multiply numbers on a computer? What if the numbers are huge and you used FFT-accelerated ...

13

This is fixed in version 9. This came up on MathGroup before. Since it hasn't been fixed for so long, I wasn't sure if it was really a bug, so I did some spelunking (and some speculation) today to find out what's happening. To jump to the end: I think it's a bug. First, let's see what arguments does LogLinearPlot really pass to the function: ...

13

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[NDSolveProcessEquations[{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 *) ...

13

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

13

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

12

If you know the equation defining your ellipsoid you could use Boole[] to constrain the integration domain : myF[x_,y_]=Abs[x+y] NIntegrate[Boole[(x/3)^2 + (y/2)^2 <= 1] myF[x,y], {x, -5, 5}, {y, -5, 5}] Note that this will actually prevent myF[x, y] from being evaluated outside the domain specified by Boole. This feature of NIntegrate is described ...

12

This approach finds equilibrium by checking that all derivatives up to the order of the differential equation are below a threshold. Following the template (defined below) suggested by the OP, here is an example for a damped harmonic oscillator: Needs["DifferentialEquationsInterpolatingFunctionAnatomy`"]; eqns1 = {a''[t] == Pi^2/2500 - (Pi^2*a[t])/2500 - ...

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