19

Your expression can be written in a more numerically stable form: Sqrt[x + 2] - Sqrt[x] == 2/(Sqrt[x + 2] + Sqrt[x]) // FullSimplify (* True *) This evaluates without significant loss of precision. 2/(Sqrt[x + 2] + Sqrt[x]) /. x -> N[32*^15] (* 5.59017*10^-9 *) 2/(Sqrt[x + 2] + Sqrt[x]) /. x -> N[32*^15, 30] (* 5.59016994374947415367652880074*10^-9 ...


18

As was pointed out above, this is a good summary of Mathematica's constrained optimization methods. Read through this if you want to know a lot more. A quick answer is below: The answer to your question is strongly dependent on the function you want to maximize. Convex functions can be maximized quite easily, with the error controlled by the PrecisionGoal ...


15

We can adapt the MonitorMethod: Options[MonitorMethod] = {Method -> Automatic, "MonitorFunction" -> Function[{h, state, meth}, Print[{"H" -> h, "SD" -> state@"SolutionData"}]]}; MonitorMethod /: NDSolve`InitializeMethod[MonitorMethod, stepmode_, sd_, rhs_, state_, OptionsPattern[MonitorMethod]] := Module[{submethod, mf}, ...


14

Can also use exact numbers. f[x_] = Cos[x] - E^(-27 x/10); f[17 10^-26]//N[#,50]& (*4.5899999999999999999999988020950000000000000000002*10^-25*) Another way is to increase the floating point digits. f[x_] = Cos[x] - E^(-2.7`50 x) f[1.7`50 10^-25] (*4.5899999999999999999999988020950000000000000000002*10^-25*) These methods pretty much work for ...


13

You can do 5.291`4/0.003`4 (* 1764. *) Precision[%] (* 3.69897 *) But as you noted, the precision of the result is lower than 4 due to precision tracking. Here's how to turn off precision tracking: How to disable roundoff error tracking in arbitrary precision arithmetic? Since you want to do this for demonstration purposes, I suggest using the ...


13

Some experimentation shows that the NDSolve problem is associated exclusively with z with no feedback to x or yT, because dz and dyT are independent of z. It is, therefore, convenient first to compute x and yT, and only then to compute z. For paramsGood, deqns = {dx, dyT, initx, bc1x, bc2x, inityT}; deqnz = {dz, initz, bc1z, bc2z}; {xSoln, yTSoln} = ...


12

By using arbitrary precision numbers instead of machine numbers: x = 3.2`100 10^30; Sqrt[x + 2] - Sqrt[x] 5.59016994374947424102293417182731712454783504367865447721289523170435 *10^-16 Or, even better, by using exact numbers and by doing the conversion afterwards (this forces Mathematica to compute all 100 leading digits): x = 32/10 10^30; N[Sqrt[x + ...


12

I'll add more details when I have my copy of Wagon's "Mathematica in Action" again, but as I mentioned in a comment, one possible way to check your solution would be to "integrate backwards", with initial conditions taken from your first call to NDSolve[]. Using Michael's example: sol1 = NDSolveValue[{x''[t] + x[t] == 0, x[0] == 1, x'[0] == 1}, x, {t, 0, 2},...


10

First convert the expression to trigonometric form: y = Cos[x] - Exp[-2.7*x] // ExpToTrig Cos[x] - Cosh[2.7 x] + Sinh[2.7 x] y /. x -> 1.7*10^-25 4.59*10^-25 This is the exact solution. Note: Building on the interesting comments on whether the above is the exact solution, it is worth noting that there is an exact relationship between the three ...


10

I'll mimic the precision/accuracy handling for FindRoot as indicated in its documentation: The default settings for AccuracyGoal and PrecisionGoal are WorkingPrecision/2. The setting for AccuracyGoal specifies the number of digits of accuracy to seek in both the value of the position of the root, and the value of the function at the root. The ...


10

To give a simple answer to the question in the title: There is dedicated hardware in your CPU for machine precision computations. In contrast, arbitrary precision computations have to be emulated in software. That is why machine precision computations are orders of magnitude faster. And that is also why Mathematica uses machine precision as default. Most of ...


10

When you use Evaluate, the PDF evaluates to an exponential: pdf = PDF[NormalDistribution[1.6, .2], x] 1.99471 E^(-12.5 (-1.6 + x)^2) When you evaluate this object with a large enough number x (on the order of 10), the resulting object is too small to represent as a machine number: pdf /. x -> 10 General::munfl: Exp[-882.] is too small to ...


10

For this problem, you must specify a solution method. Since the system of equations is nonlinear and the equation for y does not contain derivatives with respect to r, we must choose "MethodOfLines". It solves all problems right away. The code contains Quiet, since function xI is used outside the machine number, for example, 1/2.71828^726.112is calculated. ...


10

Use the asymptotic form of the function, i.e. Normal@Series[Sqrt[x + 2] - Sqrt[x], {x, ∞, 2}] % /. x -> 3.2*10^30 (* -(1/(2 x^(3/2))) + 1/Sqrt[x] *) (* 5.59017*10^-16 *) Note that we don't even need the second term to this level of precision: Normal@Series[Sqrt[x + 2] - Sqrt[x], {x, ∞, 1}] % /. x -> 3.2*10^30 (* 1/Sqrt[x] *) (* 5.59017*10^-16 *)


9

Map[Precision, tt2, {2}] shows the Gram-Schmidt process loses about 33 digits of precision on the OP's example. This precision loss estimate is the result of Mathematica estimating the propagated error and is normally an upper bound on the error. When the error is greater than the point-estimate of the result, the result is returned as an arbitrary-...


9

Because of the way you are generating the lists. What you are plotting is just a list of "y" values without the corresponding "x", that with ListPlot are plotted by Mathematica from 1 to n (with n the number of elements). If you start the list from a "x" value larger than one, Mathematica still plots from 1 to n. Try with the following code where I ...


9

Fully NDSolve-based Solution Adjustion for spatial step size together with temporal step size helps. I've used parameters mentioned in the comment for testing: mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, ...


8

At least for the OP's special case, there is a simpler way to calculate the lower triangular matrix tt: tt = LowerTriangularize[Inverse[Transpose[CholeskyDecomposition[N[mat, 20]]]]]; Norm[tt.mat.Transpose[tt] - IdentityMatrix[Length[mat]], "Frobenius"] 0 To avoid the cleanup needed to be done with LowerTriangularize[], here is an equivalent method (...


8

Actually the code is already broken in v11.3, the underlying issue is, "CatchMachineUnderflow" option is removed in this version. To fix the code, just turn to arbitrary precision and choose a high enough WorkingPrecision: SimulatePendulum[force_] := AnimatePendulum[ First[NDSolve[ SetPrecision[#, Infinity] &@{eqns /. f[t] -> force, {...


7

Because the precision necessary to represent 18154980120849865. slightly exceeds, $MachinePrecision (* 15.9546 *) RiemannSiegelZ gives an inaccurate answer, RiemannSiegelZ[18154980120849865.] (* -0.563204 *) as compared with RiemannSiegelZ[SetPrecision[18154980120849865, 30] (* 1.22954136 *) and the same is true of the entire ListLinePlot shown in the ...


7

The problem does not seem to be in Plot but in NDEigensystem. Apparently, the default method used for your function is not ideal. If you provide a method explicitly it seems to work better. {vals, funs} = NDEigensystem[ -(A/x)*D[x*D[ψ[x], {x}], {x}] + V[x]*ψ[x], ψ[x], {x, 0, 5}, num, Method -> {"PDEDiscretization" -> ...


7

One idea I like to use for these kinds of integrals is to add an auxiliary variable and a Dirac delta function, convert the Dirac delta function to it's integral formulation, and then do a bunch of simple 1D integrals. For your case, this would proceed as follows, starting from: $$\underset{t\in \mathbb{R}^n}{\int }e^{-t_1^4-t_2^4-\ldots \ -t_n^4-\left(1-...


7

How can I make it return much more precise results? There is not a great discussion of precision and accuracy control, AFAIK. Familiarize yourself with the functions in this guide: Precision & Accuracy Control Take care to distinguish distinct but related ones: N and SetPrecision MachinePrecision and $MachinePrecision Precision and Accuracy f[N[x, p]...


7

Replace your approximate input with exact rational numbers, reduce the result to a single exact algebraic number, and evaluate that numerically. f[xx_] := With[{x = Rationalize[xx, 0]}, N[RootReduce[Sqrt[x + 2] - Sqrt[x]]]] f[3.2 10^30] (* 5.59017*10^-16 *)


7

The primary issue in this question is plotting resolution. Both Plot and ContourPlot sample the function to be plotted at a finite number of locations. If none of these locations is sufficiently near a portion of the desired curve, then that portion may be missed. In this case, greatly increasing PlotPoints and setting MaxRecursion yields, Plot[function[3,...


6

Use FindRoot directly with arbitrary-precision n = 5; n1 = 4; n2 = 6; γ = 1/20; α = 1/370; icdf = x /. FindRoot[CDF[NoncentralFRatioDistribution[1, n1 - 1, n1/γ^2], x] == 1 - α - (n - n1)/n2, {x, 500}, WorkingPrecision -> $MachinePrecision] (* 5468.146403807255 *) Verifying, (CDF[NoncentralFRatioDistribution[1, n1 - 1, n1/γ^2], icdf] // ...


6

Since you liked Szabolcs' answer, but don't want to type ComputerNumber when entering a number, you can hack the input to Mathematica using $PreRead to automatically interpret every number in the input as a ComputerNumber: <<ComputerArithmetic` ComputerArithmeticOn[] := ( $PreRead = ( # /. s_String /; SyntaxQ@s && Head@ToExpression@s ===...


6

As noted in the Comments, this integration is plagued by precision problems. To proceed, factor the huge constant (h/(((1.989)*10^(30)))*(r^3)) u (* 3.44312*10^51 *) from the integrand and then FullSimplify and Rationalize the functions F[x, y]*G[x, y]] rat = Rationalize[FullSimplify[F[x, y]*G[x, y]], 0]; in order to achieve any possible cancellations ...


6

You can use Congruent. ClearAll[Congruent] TOL = 0.05; Congruent[a_, b_] := If[Abs[a - b] > 0, Abs[a - b] / Norm[{a, b}, Infinity] <= TOL, True]; Addendum your answer does not show how I test these equations with Congruent. Can you show me how I use Congruent? Please note that lines 37 and 39 on the attached screenshot use Congruent. It can ...


6

Some more ways, with the relative error e = 0.05: Block[{Internal`$EqualTolerance = MachinePrecision + Log10[e]}, {x + y == 250., z + p == 65.} (* advantage: equations written in terms of == *) ] (* {True, False} *) {SetPrecision[x + y, -Log10[e]] == SetPrecision[250, -Log10[e]], SetPrecision[z + p, -Log10[e]] == SetPrecision[65, -Log10[e]]} (* {...


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