Hot answers tagged

13

Numerics in Mathematica can be as precise as you like. However, precision comes at price; you pay for it in computation time and in additional coding effort. In Mathematica there are several computational classes of non-complex numbers, which form a tree like this. The computation you made was made with machine reals because you included 0.5 as a term. ...


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


7

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


5

The problem seems to be that you give the initial condition at a different point that your lower limit of the integration. If you change x1 by 0 things seem to work: (*INPUT*)MP = 1; m = 0.2 MP; c = Sqrt[3/2]; x1 = 0; x2 = 4000; (*ODE*) Clear[s2] Table[s2[i] = NDSolve[{y''[x] + c (Sqrt[m^2 y[x]^2 + y'[x]^2]) y'[x] + m^2 y[x] == 0, y[0] == ...


4

The main step is translating the notation of the question to a matrix problem that Mathematica can solve numerically. Once you know what the symbols in the question actually mean, the rest is quite straightforward. Here I define the matrices a and ketBra corresponding to the physical interpretation of your symbols, and then rewrite the Hamiltonian using ...


4

First of all there is somewhere in Mma a package for numerical calculation of derivatives, but I did not manage to find a reference. To offer a way to calculate the derivative. You could use the interpolation function. Here is your list: lst = {{0, 1}, {5/999, 0.999925}, {5/333, 0.9997}, {25/999, 0.999325}, {35/999, 0.998801}, {5/111, 0.998127}, ...


4

ClearAll["Global`*"] T = Table[0, {i, 1, 3}, {j, 1, 3}, {k, 1, 3}]; Part[Part[Part[T, 1], 1], 1] = E^(3 (KK + h)); Part[Part[Part[T, 1], 3], 3] = 3 E^(-KK + h); Part[Part[Part[T, 2], 2], 2] = E^(3 (KK - h)); Part[Part[Part[T, 2], 3], 3] = 3 E^-(KK + h); T = Normal[Symmetrize[T]]; h = 0; KK = 10; im = Partition[Flatten[TensorContract[TensorProduct[T, T], ...


3

If all your problems are of this sort, simply use Solve. Solve acts by manipulating the symbols, which in these cases is trivial, so that Solve[1/x == 10^(-50), x] gives the proper answer directly. Mathematica recognizes this equation as $1/x = a$ and solves, yielding $a = 1/x = 10^{(-(-50))} = 10^{50}$. In contrast, NSolve performs a numerical ...


3

Analysis of the error (bug?) We can see from the trace below that the second limit, which carries out a ratio test for the product, mistakenly yields -17 (which would indicate divergence, if correct). Trace[ NProduct[(n^2)!/stirling[n^2], {n, 1, Infinity}], _Limit, TraceInternal -> True, TraceForward -> True] There might have been some ...


3

The problem is that you are not using large enough precision goal for an oscillatory integrand with very small absolute values over the integration range. Look at the log plot -- there are more oscillations than the ones visible with Plot: The remedy is to use higher precision goal. Try this: NIntegrate[ Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/ 2 ...


2

I'd suggest a simpler approach that seems to work just fine: Use arbitrary precision input if you can. In your case, instead of d -> 0.001/5, use d -> 5/1000 integrand = Re[(x^4 Sin[x b])/(9 + 4 d^2 x^2)^6 1/2 E^(-(1/2) x^2 - I x c - 1/2 a^2) (Erfc[(I x - c)/Sqrt[2]] + E^(2 I x c) Erfc[(I x + c)/Sqrt[2]]) /. {b -> 5, ...


2

Short answer: One workaround is to use Method -> "LevinRule" instead. Long answer: As mentioned by Xavier in a comment above, changing BesselJ[0, x] to BesselJ[0, Re[x]] resolves the issue: NIntegrate[{1, -1} BesselJ[0, Re@x], {x, 0, ∞}, WorkingPrecision -> 32, Method -> "ExtrapolatingOscillatory"] Precision /@ % ...


2

Is there any reason you can't solve for logarithms? If you do, NSolve handles the provided example pretty straightforwardly: Thread[y -> (Exp[y] /. NSolve[-y == Log[10^(-50)], y])] (* {y -> 1.*10^50} *)


2

A point has measure 0 in any numeric approximation possible, so it's not only mathematica who doesn't like point-like conditions. soln = NDSolveValue[{Laplacian[u[x, y], {x, y}] == NeumannValue[-1, y == -1. && Abs[x] < 0.1] + NeumannValue[1, y == 1 && Abs[x] < 0.1], u[0, y] == 0}, u, {x, y} \[Element] Rectangle[{-1., ...


2

In version 10.2 or later you can use the EXPERIMENTAL function FindFormula lst = {{0, 1}, {5/999, 0.999925}, {5/333, 0.9997}, {25/999, 0.999325}, {35/999, 0.998801}, {5/111, 0.998127}, {55/999, 0.997305}, {65/999, 0.996335}, {25/333, 0.995217}, {85/999, 0.993952}, {95/999, 0.992542}, {35/333, 0.990986}, {115/999, 0.989287}, {125/999, ...


2

As I noted in a comment above, it is quite possible for the solution of a nonlinear PDE to become singular at finite t, and that appears to be occurring here. And, as noted by Dr. Wolfgang Hintze, the right side of the PDE, when integrated over {x, -6, 6} is zero. So, the integral of f must be a constant, and indeed it is. For the parameters given in the ...


1

This is just a long comment really. Out of curiosity, I compared the result of three packages: Mathematica, MATLAB and eig_sym from Armadillo (compiled on OS X). You said that your C++ code uses Armadillo. I get very close but not identical results: Mathematica MachinePrecision: {1.08568*10^15, -1.08568*10^15, 2.04979*10^13, 1.98037*10^13, ...


1

I think that the function f in Mr. Wizard's answer linked above may be what you want. But let me ask you, when you say you want N[1/Sqrt[3], 1] which is 0.6, do you want that precision maintained throughout calculations? That is, if I enter 1/.6 it returns 1.66667 with a repeating decimal representation. But if I force the precision to be 1 on this ...



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