# Tag Info

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

25

The default value of $NumberMarks Automatic means that  should by default be used in arbitrary-precision but not machine-precision numbers. Arbitrary-precision numbers can contain an arbitrary number of digits e.g. : Sqrt[321] == 1.73205080756887729353 Machine numbers contain the same number of digits and maintain no information on their ... 24 The backtick is a short-hand to mark the precision of your output. If it is not followed by any number, it denotes machine precision. You can denote arbitrary precision by including a number, as for example, 0.320. By default, these are not displayed in StandardForm, which is why you see them only when copying, at which point it gets converted to InputForm.... 21 I have figured out why you are getting the structure you are getting. The reason has to do with your initial choice of the angle, which you set at$0$in the Nest[] statement. The actual image is generated by choosing the mean result of iterating the map for many initial values chosen uniformly at random in$[0,1]$. With$n = 50$iterations and$m = 20$... 19 If you calculate Log[2,Log[2,$MaxNumber]], you'll get 29.999999828017338886225739 which is remarkably close to 30. Therefore I conclude that Mathematica calculates with a 31-bit exponent (1 bit for the exponent's sign). Which means that if Mathematica uses the same ordering as IEEE floats (i.e. first sign bit, then exponent, then mantissa), the first 32 ...

17

As the comments indicate, there is no completely hardware-based solution - but that doesn't mean you can't do some tweaking. The trick is always: stick with machine precision as long as you can, then switch to arbitrary precision only to refine your results. Instead of making up an example (which is hard because Mathematica implements the above principle ...

17

David Goldberg ("What every computer scientist should know about floating-point arithmetic", ACM Computing Surveys, Vol 23, No 1, March 1991, p 12, Th 4) gives pseudocode that is equivalent to log1p[x_Real] := With[{w = 1 + x}, If[w-1 == 0, x, x*Log@w/(w-1)]] EDIT - Following Mark Adler's comments, I checked the binary representation of the results (using ...

16

As it seems to depend on more than machine bits I'm curious what $MaxNumber various Mathematica installs have. If your setup is different please fill in system information and Log2 @ Log2 @$MaxNumber // Round in the table below: $$\begin{array}{r|c|c|l|c} \text{OS} & \text{Bits} & \text{Version} & \text{\MaxNumber} & \log_2\log_2\\ \... 15 I wonder whether I have understood your question correctly because I know you'll be aware of Clip data = Clip[#, {-MaxMachineNumber, MaxMachineNumber}] & /@ {0, Exp[1000.]} (* ==> {0, 1.797693135*10^308} *) Precision /@ data (* ==> {\[Infinity], MachinePrecision} *) data = RandomReal[10, {10, 2}]~Join~{{0, Exp[1000.]}}; Graphics[Point[... 14 How do I tell mathematica that all numbers e.g. 1.5 are actually 20 Digits precision? SetPrecision on all numbers or add the 20 everywhere? You could force this with PreRead. This naive definition is likely inefficient and probably breaks a number of corner cases I have not considered, but here is a rough demonstration: PreRead = (# /. ... 14 I don't think Mathematica has that function. Seems like it should. (Same for expm1().) You should not need to resort to non-machine arithmetic to get the right answer. Here is something that will do the trick using only machine arithmetic, if the input is a machine number: log1p[x_] := If[MachineNumberQ[x], If[x < 0.5, If[# - 1 == 0, x, x ... 11 As others have mentioned, the wrong result is given by N[expr] and the errors are due to cancellation. Let's discuss a bit why N[expr, 3] is able to give a good result. Mathematica can do computations with inexact ( = floating point) numbers in two ways: Using the computers native floating point arithmetic, which is very fast, but has no precision ... 11 The OP linked to an answer of mine for interpolating over general point sets; for constructing a single interpolating function, a slight modification of my procedure is needed. (In particular, you don't need centripetal or chord-length parametrization in this case.) Let's start with some data: data = {{0, 0}, {1/10, 3/10}, {1/2, 3/5}, {1, -1/5}, {2, 3}, {3,... 10 This question looks as a duplicate of these questions: How to create internally optimized expression for computing with high WorkingPrecision? How to work with ExperimentalNumericalFunction? An internally optimized version of the original function can be created as follows: n = 500; f = ExperimentalCreateNumericalFunction[{a, b}, Unevaluated[Nest[# ... 9 The exact equality comparison returns unevaluated. Root[1 + 6 #1 - 12 #1^2 - 32 #1^3 + 16 #1^4 + 32 #1^5 &, 5] == root[1] (* Out[900]= Root[1 + 6 #1 - 12 #1^2 - 32 #1^3 + 16 #1^4 + 32 #1^5 &, 5] == Cos[(2 \[Pi])/11] *) The numerical values agree to all significant digits. N[Root[1 + 6 #1 - 12 #1^2 - 32 #1^3 + 16 #1^4 + 32 #1^5 &, 5], 20] ... 9 Here is my implementation of quadric spline arbitrary precision interpolation which is defined exactly as in Interpolation with options Method->"Spline", InterpolationOrder->2. Theoretical background Quadric spline interpolation for n datapoints is defined as a Piecewise function which consists of n-2 parabolas which are spliced together in the ... 9 The last column seems in error. Here's a workaround for the sample problem, although it does not fix NDSolve`FiniteDifferenceDerivative: dx2 = dx1; dx2[[All, -1]] = -Reverse@dx1[[All, 1]] (* {-0.50000000000000000000, 0.25989153247414500869, -0.29289321881345247560, 0.36161567304292239214, -0.50000000000000000000, 0.80995720221088751026, -1.... 8 The loss of precision of this function by itself seems to be fairly modest: the big problem comes from the huge integers produced when expanding it. You can see that more clearly by changing the function to be tail recursive (so that larger values of n will be accessible without blowing the stack): ClearAll[z2]; z2[result_: 0, n_Integer, c_] := z2[result^2 +... 7 Your function isn't evaluating correctly when given inexact input: In[16]:= Table[f2b[N[E^-k]], {k, 0, 50, 10}] Out[16]= {-0.01730248257001, -0.01784636881283, -0.01785014954397, -0.01785017502377, -0.01785017519545, -0.01785017519660} If we force f2b to be evaluated with exact inputs (delaying the numericization of the result) we get the ... 7 The determinant of the matrix A in this case is about 10^282 The determinant isn't very useful, but the condition number is: You can use SingularValuesList to get the largest and the smallest singular value. If the ratio between the two is too large, the matrix is ill-conditioned. Solving an ill-conditioned linear system will still give "exact" results (... 7 This is not an answer. But I don't believe we should close this question as "easily found in the documentation". Numerics in Mathematica is an extremely complicated and mostly undocumented subject, where several mathematical concepts run up against each other in subtle and non-trivial ways. I have been thinking for some time that we ought to address this ... 6 Before people get any ideas: although we have the identity:$${}_2 F_1\left({{1,1}\atop{m}}\mid -1\right)=\frac{m-1}{2}\Phi\left(\frac12,1,m-1\right) where $\Phi(z,s,a)$ is the Lerch transcendent; or, in Mathematica notation: Hypergeometric2F1[1, 1, m, -1] == (m - 1) HurwitzLerchPhi[1/2, 1, m - 1]/2 the computation becomes even more unstable with that ...

6

My guess is that the wrong answer is the one given by N[expr] and not N[expr,3]. My mind is simple and I cannot manage big numbers, so, let's give 'em names: aN = 95881665812878; bN = 120000000000000; cN = 121576521638975; dN = 321097753837557; eN = Log[1/3 (1 + (19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/3))]; fN = 321097753837557; gN = Log[(-1 + 1/3 (...

6

Perhaps SetAttributes[log1p, NumericFunction]; N[log1p[x_?MachineNumberQ], _] := N@log1p[SetPrecision[x, $MachinePrecision]]; N[log1p[x_?NumericQ], {MachinePrecision, MachinePrecision}] := N@N[log1p[x],$MachinePrecision]; N[log1p[x_?NumericQ], a_] := N[Log[1 + x], a] EDIT This probably makes more sense ClearAll[log1p] log1p[x_?MachineNumberQ] := ...

6

Problem The problem with Log[1. + 1.*^-15] not yielding 1.*^-15 is not due to Log, but to MachinePrecision inputs, which I think the OP implied in the question statement: 1 + 1.*^-15 % - 1 (* 1. 1.11022*10^-15 *) So Log[1 + 1.*^-15] does return the right answer, 1.11022*10^-15, for the actual input. Solution Here is a simple way to get log1p-type ...

6

The following is a shameful plug of J. M.'s answer that you already linked to show how to get the parametrization by using BSplineBasis[] with arbitrary precision and packed into a function: splineInterp[data_, order_, prec_] := Module[{parametrizeCurve, tvals, bas, ctrlpts, knots}, parametrizeCurve[pts_ /; MatrixQ[pts, NumericQ], a : (_?NumericQ) : ...

6

For large arguments, is it not better to use Erfc? For example: myErf[z_, prec_:$MachinePrecision] := 1 - N[Erfc[z], prec]; myErf[200] (* -> 0.9999999999999999...9999999999999531064070401639850196 *) This also has the advantage of not requiring involved calculations at very high precision, allowing the answer to be produced almost immediately. 6 not an answer, just a cool animation: Looping over the initial seed values from 0-1 6 Incidentally, here is a higher resolution picture of the image produced by heropup's answer, plotted over both positive and negative values of the two circle map parameters: Blue represents negative, yellow represents positive. You can view a far larger 6001x6001 pixel image (about 29MB) here. Here is the code used to generate the upper right quarter of ... 5 Straight from the Mathematica documentation of SetPrecision bit = Log[10., 2.]; f[x_] := Module[{p = Precision[x], lx}, lx = Block[{$MaxPrecision = p, \$MinPrecision = p},4*x*(1 - x)]; SetPrecision[lx, p - bit]] then testing x0 = N[1/3, 20]; fl = NestList[f, x0, 20] As mentioned in the comment of your question, those related questions has got tons of ...

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