# Tag Info

## Hot answers tagged numerics

94

I think the reason is to prevent biasing numbers on average upward or downward. For example if you have a list of numbers that include a lot of x.5 and you were to round all these upward, the average magnitude of the list would also shift upward. Likewise if you round downward, downward. By rounding to the nearest even number these would balance out, ...

83

OK, there is good news and there is bad news. In the current version 10 there is no way to do this directly. That's the bad news. The good news is that finite element framework used within NDSolve is exposed and documented; for maximum "hackability" convenience. Let's start with a region that @MarkMcClure would consider interesting. We load our favorite ...

67

Control the Precision and Accuracy of Numerical Results This is an excellent question. Of course everyone could claim highest accuracy for their 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 ...

60

I've encapsulated the code of the mysterious user21 into a helmholzSolve command. The code is at the end of this post. It adds very little to user21's code but it does allow us to examine multiple examples quite easily, though it has certainly not been tested extensively and could be improved quite a lot I'm sure. It should be called as follows: {ev,if,...

57

It looks like you need some basic facts on numerical Fourier transforms. I am going to assume that you are going from the time domain to the frequency domain. The number of points in the time domain equals the number of points in the frequency domain. As the data is sampled in the time domain i.e. is a set of equally spaced points, then in the frequency ...

55

If you Rationalize your real numbers you will be able to use Mathematica's arbitrary precision engine: poly2 = Rationalize[poly[z], 0]; Plot[poly2, {z, 0, 1}, WorkingPrecision -> 50] Arbitrary and machine precision Mathematica has two kinds of numeric calculations: machine precision, and arbitrary precision. Machine precision is fast but is limited ...

52

Here is my latest code for this function, from Chapter 12 of the third edition of "Mathematica in Action". It is pretty short, but I will let you work out if it is faster or more robust than yours. Note the PlotPoints option for difficult cases. FindRoots2D::usage = "FindRoots2D[funcs,{x,a,b},{y,c,d}] finds all nontangential solutions to {f=0, g=0} in ...

52

In general, a typical root of a negative number is complex, so you need to get rid of most roots. A nice approach would be Root, e.g. Root[ x^3 + 8, #] & /@ Range[3] {-2, 1 - I Sqrt[3], 1 + I Sqrt[3]} To get only real roots you can do : Select[Root[ x^3 + 8, #] & /@ Range[3], Re[#] == # &] {-2} This is a handy approach when you have ...

49

It is called bankers' rounding. The rationale is that the rounding behaves "nicely" even if you have negative numbers, i.e. rounding commutes with negation, which is nice.

49

The first example seems to intentionally set Mathematica up to "fail" by specifying insufficient input accuracy. With additional precision: ClearAll[s] s[i_] := s[i] = 2*s[i - 1] - 3*s[i - 1]^2 s[0] = 0.330; s[40] 0.333333 And Mathematica is capable of far greater precision if necessary: ClearAll[s] $RecursionLimit = ∞ s[i_] := s[i] = 2*s[i - 1] - 3*s[... 43 I was recently reminded that the following "functions" are settable, and I was surprised (even though I've seen this before). So I thought that it is valuable to share this information. There are three critical properties of symbolic constants: The Constant attribute, which is used by Dt The ability to tell that it is a numeric expression (NumericQ[Pi] ===... 42 Here's some code that I used recently, based on code by Paul Abbott [1, 2]. Clear[TranscendentalRecognize] TranscendentalRecognize[num_Real, basis_List, ord_?Positive, debug:(True|False):False] := Module[{vect, mat, lr, ans}, vect = Round[10^Floor[ord - 1] Join[{num}, N[basis, ord]]]; mat = Append[IdentityMatrix[Length[vect]], vect]; lr = ... 39 Borrowing almost verbatim from a recent response about finding extrema, here is a method that is useful when your function is differentiable and hence can be "tracked" by NDSolve. f[x_] := BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]] In[191]:= zeros = Reap[soln = y[x] /. First[ NDSolve[{y'[x] == Evaluate[D[f[x], x]], y[10] == (f[10])}, ... 38 This is ContourPlot based but seems much shorter: FindCrossings2D[{f_, g_}, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}] := {x, y} /. (FindRoot[{f[x, y] == 0, g[x, y] == 0}, {{x, #[[1]]}, {y, #[[2]]}}] & /@ (ContourPlot[{f[x, y] == 0, g[x, y] == 0}, {x, xmin, xmax}, {y, ymin, ymax}][[1, 1]])) It works: f[x_, y_] := -Cos[y] + 2 y Cos[y^2] Cos[2 x]; ... 38 Oleksandr is correct about the way evaluation works. a/b seems to be interpreted (parsed) directly as Times[a, Power[b,-1]], or more readably:$a\times b^{-1}$. Divide[a,b] is interpreted as is. Evaluation then proceeds from these forms, and the arithmetic is carried out differently for the two cases: either$a\times (1/b)$or$a/b$. Here are some ... 38 Without commenting on how much attention one should pay to marketing literature: the second example is somewhat relevant. Mathematica's polynomial factoring algorithm is known to be at least fifteen years behind the state of the art, and things that Maple will factor in seconds will go away (literally) forever in Mathematica. This is, of course, not too ... 36 I can offer a round-about method. First compute the numerical approximation. I obtain, to high precision, In[24]:= N[Sum[1/(2*n!), {n, 0, 100}], 100] Out[24]= 1.\ 3591409142295226176801437356763312488786235468499797874834838138620383\ 15176773797285691089262583214 Now paste that into a Wolfram|Alpha query, accessed by clicking on the '+' sign at upper ... 35 You have selected the wrong digit. Mathematica gets the digit in the million-th decimal place right if the calculation is performed correctly. q = N[Pi, 1000010]; RealDigits[q][[1, 1000001]] 1 I take the 1000001-th digit because RealDigits includes the integer part, 3. Update It is really important to use RealDigits to decide this question. Looking at ... 34 Let me give a different approach. FindRoot does a good job, but maybe we can calculate the seed-points in a different way. When you want to find the common roots of$f(x,y)$and$g(x,y)$you can transform the problem into one equation which has the same roots $$0=f(x,y)^2+g(x,y)^2$$ The nice property here, which I will use is that the right hand side of this ... 34 Unfortunately, this requires a lot more devious trickery than I would have preferred. As noted in the documentation at tutorial/UnconstrainedOptimizationQuasiNewtonMethods, the Hessian is not formed directly in the BFGS method, so we have to recover it from the Cholesky factors. However, all of this is done inside the kernel where we cannot access it using ... 34 Obviously, for large negative inputs, Exp will produce very small numbers. While this isn't intrinsically problematic, it so happens that, by default, Mathematica deals with machine underflow by converting the affected values to an arbitrary precision representation in order to avoid catastrophic loss of precision. However, sometimes one would rather ... 34 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 ... 33 The critical issue in the first example is that Mathematica is using significance arithmetic to track precision. This is certainly billed as a feature by Wolfram Research. As we see in this example though, it can be portrayed as a weakness. In truth, you might need to know what you're doing to use it correctly. In this answer, I mentioned that significance ... 32 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.... 32 If this is something you want in general, try: SetOptions[$FrontEnd, PrintPrecision-> 10] and if you just want it for a specific notebook, then do: SetOptions[InputNotebook[], PrintPrecision-> 10]

32

To create ExperimentalNumericalFunction, one needs to evaluate ExperimentalCreateNumericalFunction[vars, expr, dims] where vars is a list of arguments, expr - the expression from which the numerical function will be created, dims - the dimensions of the output matrix produced by this expression. If the output is scalar, then dims should be set to {}. It ...

31

There are two rather different scenarios for numerical derivatives: Differentiating a continuous function that's only defined numerically Approximating the derivative of a list of data that could itself be generated numerically For scenario 1, here is an example function and its derivative: f[x_?NumericQ] := BesselJ[1, x] Needs["NumericalCalculus"] ...

31

Following the advice in comments, I've made a test library for BesselJ[1, #] & function to evaluate via GSL. I still consider it a workaround, so if you find a way to use Mathematica built-in functions with good performance, please do make a new answer. Needs["CCompilerDriver"] besselJ1src = " #include \"WolframLibrary.h\" DLLEXPORT mint ...

31

Introduction The first section below can be found in standard numerical analysis textbooks. Most current textbooks seem to assume a working environment such as MATLAB or a programming language such C, Python, etc. in which there are at most two choices for working precision: single and double. As a result, WorkingPrecision as a parameter is not much ...

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