Working with Precision in Mathematica

Jayanta Phadikar, Wolfram Technology Engineer

Join the Conversation #WolframTechConf

High-precision arithmetic

Traditional programming languages incorporate limited-precision arithmetic model

Examples: Single precision (32 bits), Double precision (64 bits)

For single precision, computer can use maximum 32 bits to represent a number

Arbitrarily large information about a number can’t be stored e.g.

136891479058588375991326027382088315966463695625337436471480190078368997177499076593800206155688941388250484440597994042813512732765695774566001

5.58620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862

High precision arithmetic is useful when:

i. The result of a computation is highly sensitive to the initial conditions or intermediate results

ii. We desire very accurate numerical results

Mathematica allows arbitrary-precision computation.

Computing square root of 2 up to 100 digits

N[Sqrt[2],100]

Computing Sin[10^40]

The above result is completely off.

Correct result using Mathematica ArbitraryPrecision computation:

N[Sin[10^40],10]

MachinePrecision, ArbitraryPrecision and InfinitePrecision

Three types of computation systems:

1. MachinePrecision - less accurate and faster
Error is not tracked during computation

2. ArbitraryPrecision - more accurate and slower
Error is tracked

3. InfinitePrecision - exact computation
(no errors)

A = f(F, M, α, q, l, a, EJ, b, k)

A [F_,M_,α_,q_,l_,a_,EJ_,b_,k_]:=Module[{},4 EJ k (α^3) (-((Cos[(a+l) α])^2) ((Cosh[l α])^2) ((Cosh[(a+l) α])^2) ((Sin[l α])^2)+2 Cos[l α] Cos[(a+l) α] ((Cosh[l α])^2) ((Cosh[(a+l) α])^2) Sin[l α] Sin[(a+l) α]-((Cos[l α])^2) ((Cosh[l α])^2) ((Cosh[(a+l) α])^2) ((Sin[(a+l) α])^2)+((Cos[l α])^2) ((Cos[(a+l) α])^2) ((Cosh[(a+l) α])^2) ((Sinh[l α])^2)+2 ((Cos[(a+l) α])^2) ((Cosh[(a+l) α])^2) ((Sin[l α])^2) ((Sinh[l α])^2)-2 Cos[l α] Cos[(a+l) α] ((Cosh[(a+l) α])^2) Sin[l α] Sin[(a+l) α] ((Sinh[l α])^2)+2 ((Cos[l α])^2) ((Cosh[(a+l) α])^2) ((Sin[(a+l) α])^2) ((Sinh[l α])^2)+((Cosh[(a+l) α])^2) ((Sin[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[l α])^2)-2 ((Cos[l α])^2) ((Cos[(a+l) α])^2) Cosh[l α] Cosh[(a+l) α] Sinh[l α] Sinh[(a+l) α]-2 ((Cos[(a+l) α])^2) Cosh[l α] Cosh[(a+l) α] ((Sin[l α])^2) Sinh[l α] Sinh[(a+l) α]-2 ((Cos[l α])^2) Cosh[l α] Cosh[(a+l) α] ((Sin[(a+l) α])^2) Sinh[l α] Sinh[(a+l) α]-2 Cosh[l α] Cosh[(a+l) α] ((Sin[l α])^2) ((Sin[(a+l) α])^2) Sinh[l α] Sinh[(a+l) α]+((Cos[l α])^2) ((Cos[(a+l) α])^2) ((Cosh[l α])^2) ((Sinh[(a+l) α])^2)+2 ((Cos[(a+l) α])^2) ((Cosh[l α])^2) ((Sin[l α])^2) ((Sinh[(a+l) α])^2)-2 Cos[l α] Cos[(a+l) α] ((Cosh[l α])^2) Sin[l α] Sin[(a+l) α] ((Sinh[(a+l) α])^2)+2 ((Cos[l α])^2) ((Cosh[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[(a+l) α])^2)+((Cosh[l α])^2) ((Sin[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[(a+l) α])^2)-((Cos[(a+l) α])^2) ((Sin[l α])^2) ((Sinh[l α])^2) ((Sinh[(a+l) α])^2)+2 Cos[l α] Cos[(a+l) α] Sin[l α] Sin[(a+l) α] ((Sinh[l α])^2) ((Sinh[(a+l) α])^2)-((Cos[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[l α])^2) ((Sinh[(a+l) α])^2))]

MachinePrecision

A[1000.0,0.0,((3/500)^(1/4)),0.0,50.0,50.0,50000000.0,40.0,1200000.0]

A[1000.0,0.0,((3/500)^(1/4)),0.0,50.0,50.0,50000000.0,40.0,1200000.0]//RepeatedTiming//First

ArbitraryPrecision

N[#,10]&/@A[1000,0,((3/500)^(1/4)),0,50,50,50000000,40,1200000]

N[#,10]&/@A[1000,0,((3/500)^(1/4)),0,50,50,50000000,40,1200000]//RepeatedTiming//First

InfinitePrecision

A[1000,0,((3/500)^(1/4)),0,50,50,50000000,40,1200000]

Concepts in Statistics

“Accuracy” and “Precision” concepts in statistics are different from those in Mathematica

Precision: random error
Accuracy: systematic error

https://en.wikipedia.org/wiki/Accuracy_and_precision

Enforcing Machine/Arbitrary/Infinite Precision

To enforce (for numeric functions):

MachinePrecision framework - feed the inputs in MachinePrecision

ArbitraryPrecision framework - feed the inputs in ArbitraryPrecision/”specified precision”

InfinitePrecision framework - feed exact inputs

x Precision in → x Precision out

For mixed inputs:

If there is at least one MachinePrecision input, MachinePrecision framework will be selected.
If there is at least one ArbitraryPrecision input and no MachinePrecision inputs, ArbitraryPrecision framework will be selected.

MachinePrecision

MachinePrecision numbers: 17 digits or less (without backtick)

2.3//Precision

1.23456789//Precision (* 9 digits *)

{1.2345678901234567//Precision (* 17 digits *),1.2345678901234567890123456789012//Precision (*32 digits*)}

1.2345678901234567890123456789012`//Precision

SetPrecision[1.2345678901234567890123456789012,MachinePrecision]//Precision

MachinePrecision In → MachinePrecision Out

{5.0+2.3,
6.7*(1.0+2.3/4.5),
Sin[2.3],
EllipticE[0.4]}

Precision/@%

ArbitraryPrecision

ArbitraryPrecision numbers (numbers with specified precision):

x1=2.3`20;

x2=2.30000000000000000000 ;(* 19 zeros, 21 digits *)

x3=SetPrecision[2.3,20];

Precision/@{x1,x2,x3}

ArbitraryPrecision In → ArbitraryPrecision Out

{5`20+23`20/10`20,
6.7`20*(1.0`20+2.3`20/4.5`20),
Sin[23`20/10`20],
EllipticE[0.4`20]}

Precision/@%

InfinitePrecision

InfinitePrecision numbers:

x1=23/10;

x2=45;

x3=Pi;

Precision/@{x1,x2,x3}

InfinitePrecision In → InfinitePrecision Out

{5+23/10,
67/10*(1+23/45),
Sin[23/10],
EllipticE[2/5]}

Precision/@%

Numeric expressions

Numeric expressions return a number as output
x Precision in → x Precision out

ref/NumericQ

List of all numeric functions

Row[Column[#]&/@Partition[Select[Names["System`*"],MemberQ[Attributes[#],NumericFunction]&],70],"   "]

Examples

{CubeRoot[2.0],LogGamma[2.5],Sinh[2.5],Gamma[2.5],SphericalBesselJ[2,3.5]}

Precision/@%

{CubeRoot[2.0`10],LogGamma[2.5`10],Sinh[2.5`10],Gamma[2.5`10],SphericalBesselJ[2,3.5`10]}

Precision/@%

{CubeRoot[2],LogGamma[5/2],Sinh[5/2],Gamma[5/2],SphericalBesselJ[2,7/2]}

Precision/@%

Obtaining numerical value for an expression using N[]

InfinitePrecision computation results in exact expressions. N[] can be used to determine the numerical values

N[expr] gives expr in MachinePrecision

N[expr, n] attempts to give expr with n-digit precision - ArbitraryPrecision

x=1+2/3+5/7

{N[x],N[x]//Precision}

{N[x,20],N[x,20]//Precision}

{CubeRoot[2],LogGamma[5/2],Sinh[5/2],Gamma[5/2],SphericalBesselJ[2,7/2]}

{N[CubeRoot[2]],N[LogGamma[5/2]],N[Sinh[5/2]],N[Gamma[5/2]],N[SphericalBesselJ[2,7/2]]}

Sin[7]+Cos[10]+Tan[70]

N[%]

A [F_,M_,α_,q_,l_,a_,EJ_,b_,k_]:=Module[{},4 EJ k (α^3) (-((Cos[(a+l) α])^2) ((Cosh[l α])^2) ((Cosh[(a+l) α])^2) ((Sin[l α])^2)+2 Cos[l α] Cos[(a+l) α] ((Cosh[l α])^2) ((Cosh[(a+l) α])^2) Sin[l α] Sin[(a+l) α]-((Cos[l α])^2) ((Cosh[l α])^2) ((Cosh[(a+l) α])^2) ((Sin[(a+l) α])^2)+((Cos[l α])^2) ((Cos[(a+l) α])^2) ((Cosh[(a+l) α])^2) ((Sinh[l α])^2)+2 ((Cos[(a+l) α])^2) ((Cosh[(a+l) α])^2) ((Sin[l α])^2) ((Sinh[l α])^2)-2 Cos[l α] Cos[(a+l) α] ((Cosh[(a+l) α])^2) Sin[l α] Sin[(a+l) α] ((Sinh[l α])^2)+2 ((Cos[l α])^2) ((Cosh[(a+l) α])^2) ((Sin[(a+l) α])^2) ((Sinh[l α])^2)+((Cosh[(a+l) α])^2) ((Sin[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[l α])^2)-2 ((Cos[l α])^2) ((Cos[(a+l) α])^2) Cosh[l α] Cosh[(a+l) α] Sinh[l α] Sinh[(a+l) α]-2 ((Cos[(a+l) α])^2) Cosh[l α] Cosh[(a+l) α] ((Sin[l α])^2) Sinh[l α] Sinh[(a+l) α]-2 ((Cos[l α])^2) Cosh[l α] Cosh[(a+l) α] ((Sin[(a+l) α])^2) Sinh[l α] Sinh[(a+l) α]-2 Cosh[l α] Cosh[(a+l) α] ((Sin[l α])^2) ((Sin[(a+l) α])^2) Sinh[l α] Sinh[(a+l) α]+((Cos[l α])^2) ((Cos[(a+l) α])^2) ((Cosh[l α])^2) ((Sinh[(a+l) α])^2)+2 ((Cos[(a+l) α])^2) ((Cosh[l α])^2) ((Sin[l α])^2) ((Sinh[(a+l) α])^2)-2 Cos[l α] Cos[(a+l) α] ((Cosh[l α])^2) Sin[l α] Sin[(a+l) α] ((Sinh[(a+l) α])^2)+2 ((Cos[l α])^2) ((Cosh[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[(a+l) α])^2)+((Cosh[l α])^2) ((Sin[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[(a+l) α])^2)-((Cos[(a+l) α])^2) ((Sin[l α])^2) ((Sinh[l α])^2) ((Sinh[(a+l) α])^2)+2 Cos[l α] Cos[(a+l) α] Sin[l α] Sin[(a+l) α] ((Sinh[l α])^2) ((Sinh[(a+l) α])^2)-((Cos[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[l α])^2) ((Sinh[(a+l) α])^2))]

MachinePrecision:

A[1000.0,0.0,((3/500)^(1/4)),0.0,50.0,50.0,50000000.0,40.0,1200000.0]

ArbitraryPrecision:

N[#,10]&/@A[1000,0,((3/500)^(1/4)),0,50,50,50000000,40,1200000]

Cautions and recommendations

Check with ArbitraryPrecision while working with very large/small numbers

MachinePrecision results can be completely wrong

We can always trust ArbitraryPrecision results

Example 1

MachinePrecision:

2-3*(5.2917721067*10^-9)/(.9*10^-13)+3*((5.2917721067*10^-9)/(.9*10^-13))^3-3*(5.2917721067*10^-9)*((.9*10^-13)+(5.2917721067*10^-9))^2*Exp[-2 (.9*10^-13)/(5.2917721067*10^-9)]/(.9*10^-13)^3

ArbitraryPrecision:

x=2-3*((Rationalize[5.2917721067*10^-9,0])/(Rationalize[.9*10^-13,0]))+3*((Rationalize[5.2917721067*10^-9,0])/(Rationalize[.9*10^-13,0]))^3-3*Rationalize[5.2917721067*10^-9,0]*(Rationalize[.9*10^-13,0]+Rationalize[5.2917721067*10^-9,0])^2*Exp[-2 (Rationalize[.9*10^-13,0])/(Rationalize[5.2917721067*10^-9,0])]/(Rationalize[.9*10^-13,0])^3
N[x,10]

Example 2

expr=95881665812878-120000000000000 (121576521638975-(321097753837557 Log[1/3 (1+(19-3 Sqrt[33])^(1/3)+(19+3 Sqrt[33])^(1/3))])/Log[5])+(321097753837557 Log[(-1+1/3 (1+(19-3 Sqrt[33])^(1/3)+(19+3 Sqrt[33])^(1/3)))/(-6+4/3 (1+(19-3 Sqrt[33])^(1/3)+(19+3 Sqrt[33])^(1/3)))])/Log[5];

MachinePrecision:

N[expr]

ArbitraryPrecision:

N[expr,10]

Example 3

f[n_]:=10^60*(n+10^-57-1)

MachinePrecision:

f[1.0]

ArbitraryPrecision:

f[1]

As far as possible, use InfinitePrecision computation followed by an ArbitraryPrecision computation

Practical advantages:

i. A pre InfinitePrecision computation avoids the precision loss issue

ii. We obtain a result with desired precision

iii. The procedure is simple to implement

Example 1

N[Sin[10^40.0]] (* MachinePrecision; wrong result *)

Sin[10^40`50] (* ArbitraryPrecision; precision loss: 50 -> 8 *) 
Precision[%]

Sin[10^50`50] (* ArbitraryPrecision; complete precision loss *)

N[Sin[10^50],50] (* InfinitePrecision followed by ArbitraryPrecision; we get what we request *)
Precision[%]

Example 2

(* The Rationalize[] function *)
{Rationalize[7.9],Rationalize[2.23482748392743,0]}
N[%]

2-3*(5.2917721067*10^-9)/(.9*10^-13)+3*((5.2917721067*10^-9)/(.9*10^-13))^3-3*(5.2917721067*10^-9)*((.9*10^-13)+(5.2917721067*10^-9))^2*Exp[-2 (.9*10^-13)/(5.2917721067*10^-9)]/(.9*10^-13)^3
(* MachinePrecision; wrong result *)

x=2-3*((Rationalize[5.2917721067*10^-9,0])/(Rationalize[.9*10^-13,0]))+3*((Rationalize[5.2917721067*10^-9,0])/(Rationalize[.9*10^-13,0]))^3-3*Rationalize[5.2917721067*10^-9,0]*(Rationalize[.9*10^-13,0]+Rationalize[5.2917721067*10^-9,0])^2*Exp[-2 (Rationalize[.9*10^-13,0])/(Rationalize[5.2917721067*10^-9,0])]/(Rationalize[.9*10^-13,0])^3
N[x,20]
Precision[%]
(* InfinitePrecision followed by ArbitraryPrecision; we get what we request *)

Example 3

A [F_,M_,α_,q_,l_,a_,EJ_,b_,k_]:=Module[{},4 EJ k (α^3) (-((Cos[(a+l) α])^2) ((Cosh[l α])^2) ((Cosh[(a+l) α])^2) ((Sin[l α])^2)+2 Cos[l α] Cos[(a+l) α] ((Cosh[l α])^2) ((Cosh[(a+l) α])^2) Sin[l α] Sin[(a+l) α]-((Cos[l α])^2) ((Cosh[l α])^2) ((Cosh[(a+l) α])^2) ((Sin[(a+l) α])^2)+((Cos[l α])^2) ((Cos[(a+l) α])^2) ((Cosh[(a+l) α])^2) ((Sinh[l α])^2)+2 ((Cos[(a+l) α])^2) ((Cosh[(a+l) α])^2) ((Sin[l α])^2) ((Sinh[l α])^2)-2 Cos[l α] Cos[(a+l) α] ((Cosh[(a+l) α])^2) Sin[l α] Sin[(a+l) α] ((Sinh[l α])^2)+2 ((Cos[l α])^2) ((Cosh[(a+l) α])^2) ((Sin[(a+l) α])^2) ((Sinh[l α])^2)+((Cosh[(a+l) α])^2) ((Sin[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[l α])^2)-2 ((Cos[l α])^2) ((Cos[(a+l) α])^2) Cosh[l α] Cosh[(a+l) α] Sinh[l α] Sinh[(a+l) α]-2 ((Cos[(a+l) α])^2) Cosh[l α] Cosh[(a+l) α] ((Sin[l α])^2) Sinh[l α] Sinh[(a+l) α]-2 ((Cos[l α])^2) Cosh[l α] Cosh[(a+l) α] ((Sin[(a+l) α])^2) Sinh[l α] Sinh[(a+l) α]-2 Cosh[l α] Cosh[(a+l) α] ((Sin[l α])^2) ((Sin[(a+l) α])^2) Sinh[l α] Sinh[(a+l) α]+((Cos[l α])^2) ((Cos[(a+l) α])^2) ((Cosh[l α])^2) ((Sinh[(a+l) α])^2)+2 ((Cos[(a+l) α])^2) ((Cosh[l α])^2) ((Sin[l α])^2) ((Sinh[(a+l) α])^2)-2 Cos[l α] Cos[(a+l) α] ((Cosh[l α])^2) Sin[l α] Sin[(a+l) α] ((Sinh[(a+l) α])^2)+2 ((Cos[l α])^2) ((Cosh[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[(a+l) α])^2)+((Cosh[l α])^2) ((Sin[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[(a+l) α])^2)-((Cos[(a+l) α])^2) ((Sin[l α])^2) ((Sinh[l α])^2) ((Sinh[(a+l) α])^2)+2 Cos[l α] Cos[(a+l) α] Sin[l α] Sin[(a+l) α] ((Sinh[l α])^2) ((Sinh[(a+l) α])^2)-((Cos[l α])^2) ((Sin[(a+l) α])^2) ((Sinh[l α])^2) ((Sinh[(a+l) α])^2))]

MachinePrecision:

A[1000.0,0.0,((3/500)^(1/4)),0.0,50.0,50.0,50000000.0,40.0,1200000.0]

InfinitePrecision followed by ArbitraryPrecision - much simpler

N[#,10]&/@A[1000,0,((3/500)^(1/4)),0,50,50,50000000,40,1200000]

Batch rationalization

$Pre=Function[{in},Unevaluated[in]/.r_Real?MachineNumberQ:>Condition@Rationalize[r,0],HoldAllComplete];

{5.6,3.29347829387}

3.147823*(5.8672348267*Exp[-2 (1.234672873284+6.3457)]-9.78789)
N[%]

$Pre=.

Increase $MaxExtraPrecision for N[] when needed

$MaxExtraPrecision: Maximum precision used by N[] internally
Default value: 50

N[Sin[10^100],50]

Block[{$MaxExtraPrecision=200},N[Sin[10^100],30]]

Setting up ArbitraryPrecision computation in functional operators

Default: MachinePrecision

Numerical operators on mathematical functions:

NIntegrate[], NSum[], NProduct[], FindRoot[], NSolve[], NMinimize[], FindMinimum[], NMaximize[], FindMaximum[], NDSolve[]

(tutorial/NumericalOperationsOnFunctionsOverview)

For most functions WorkingPrecision is MachinePrecision by default

s=FindRoot[x^2==2,{x,1}]

Precision[x/.s]

Options[FindRoot]

{AccuracyGoal->Automatic,Compiled->Automatic,DampingFactor->1,Evaluated->True,EvaluationMonitor->None,Jacobian->Automatic,MaxIterations->100,Method->Automatic,PrecisionGoal->Automatic,StepMonitor->None,WorkingPrecision->MachinePrecision}

All operators except NSolve[]:

(Last[Options[#]]&/@{NIntegrate, NSum, NProduct, FindRoot, NSolve, NMinimize, FindMinimum, NMaximize, FindMaximum, NDSolve})//Column

NSolve[x^2==2`50,x]

Enforcing ArbitraryPrecision

To enforce ArbitraryPrecision

i. Use all arbitrary precision inputs

ii. Set a WorkingPrecision

s1=FindRoot[x^2==2,{x,1},WorkingPrecision->50] (* Alternately PrecisionGoal *)

Precision[x/.s1]

FindMinimum[x Cos[x],{x,2},WorkingPrecision->20]

NDSolve[] Example 1

MachinePrecision - incorrect
Field-lines should go inward

deq1=x'[t]==-x[t]-0.8* y[t]+0.4 *y[t]^2;
deq2=y'[t]==2* x[t]+1.0*y[t]-x[t] *y[t]-0.2835* y[t]^2;

sol1=NDSolve[ {deq1,deq2,x[0.0]==-0.000001, y[0.0] ==0.0}, {x, y}, {t, 0.0, 5000}];

Show[ParametricPlot[ Evaluate[{x[t], y[t]} /. sol1], {t, 0.0,500},PlotRange->{{-3*10^(-6),3*10^(-6)},{-4.5*10^(-6),4.5*10^(-6)}},ImageSize->400,PlotStyle->Green],ParametricPlot[ Evaluate[{x[t], y[t]} /. sol1], {t, 0.0,20},PlotStyle->Blue]]

ArbitraryPrecision - correct

deq1R=x'[t]==-x[t]-Rationalize[0.8]* y[t]+Rationalize[0.4] *y[t]^2;
deq2R=y'[t]==2* x[t]+1*y[t]-x[t] *y[t]-Rationalize[0.2835]* y[t]^2;

sol2=NDSolve[ {deq1R,deq2R,x[0]==-Rationalize[0.000001], y[0] ==0}, {x, y}, {t, 0, 5000} ,WorkingPrecision->50,MaxSteps->200000];

Show[ParametricPlot[ Evaluate[{x[t], y[t]} /. sol2], {t, 0.0,5000},PlotRange->{{-3*10^(-6),3*10^(-6)},{-4.5*10^(-6),4.5*10^(-6)}},ImageSize->400,PlotStyle->Green],
ParametricPlot[ Evaluate[{x[t], y[t]} /. sol2], {t, 0.0,500},PlotStyle->Blue]]

Ultra-sensitive systems

In nature, we find several examples of ultra-sensitive systems - systems for which the final output is inordinately sensitive to the initial conditions

NDSolve[] Example 2

Forced vibration of a nonlinear spring-mass system with a damper

T=1000;
sol1=NDSolve[{x''[t]+15/100 x'[t]-x[t]+x[t]^3==3/10 Cos[t],x[0]==-1,x'[0]==1},x,{t,0,T},WorkingPrecision->20,MaxSteps->100000];

sol2=NDSolve[{x''[t]+15/100 x'[t]-x[t]+x[t]^3==3/10 Cos[t],x[0]==-1,x'[0]==Rationalize[1.00001,0]},x,{t,0,T},WorkingPrecision->20,MaxSteps->100000];

sol3=NDSolve[{x''[t]+0.15 x'[t]-x[t]+x[t]^3==0.3 Cos[t],x[0]==-1,x'[0]==1.00001},x,{t,0,T},MaxSteps->100000];

Time window
 t: 0 - 100

Plot[{Evaluate[x[t]/.sol1],Evaluate[x[t]/.sol2],Evaluate[x[t]/.sol3]},{t,0,T/10},PlotStyle->{Green,Blue,Directive[Black,Dashed]},ImageSize->600,PlotLegends->Placed[{"x'(0)=1; ArbitraryP","x'(0)=1.00001; ArbitraryP","x'(0}=1.00001; MachineP"},Above]]

Time window:
 t: 0 - 1000

Manipulate[Column[{Row[{Style["Time window: ",18],Style[{(Subscript[t, 0]),(Subscript[t, 0])+30},18]}],"","",Plot[{Evaluate[x[t]/.sol1],Evaluate[x[t]/.sol2],Evaluate[x[t]/.sol3]},{t,0,T},PlotStyle->{Green,Blue,Directive[Black,Dashed]},ImageSize->500,PlotRange->{{(Subscript[t, 0]),(Subscript[t, 0])+30},Automatic},PlotLegends->Placed[{"x'(0)=1; ArbitraryP","x'(0)=1.00001; ArbitraryP","x'(0}=1.00001; MachineP"},Above]]}],{(Subscript[t, 0]),0,800},SaveDefinitions->True]

Changing the display format for MachinePrecision numbers

MachinePrecision: only first 6 digits are displayed

x=N[Sqrt[10]]

N[Sqrt[10],25]

How do we display more digits in MachinePrecision?

Method 1

Style[x,PrintPrecision->15]

Method 2

NumberForm[x,15]

Method 3

InputForm[x]

Method 4

Copying and pasting

3.1622776601683795`

Largest and smallest numbers

$MaxNumber and $MinNumber

Computations with numbers greater than $MaxNumber or less than $MinNumber will fail; known as Underflow[] and Overflow[] respectively

$MaxNumber and $MinNumber are ArbitraryPrecision numbers

$MaxNumber

2 $MaxNumber

$MinNumber

$MinNumber/2

Precision[#]&/@{$MaxNumber,$MinNumber}

$MaxMachineNumber and $MinMachineNumber

Largest and smallest representable numbers in MachinePrecision:

{$MaxMachineNumber,$MinMachineNumber}

Precision[#]&/@{$MaxMachineNumber,$MinMachineNumber}

Overflow[] trapping:
During a MachinePrecision computation, if a numbers goes above $MaxMachineNumber, Mathematica will automatically switch to ArbitraryPrecision computation.

Exp[1000.0]

Precision[%]

Change in V11.3:
Mathematica doesn’t do Underflow[] trapping anymore
Developers decided to remove the feature for the sake of overall improvement in computations.

V11.2:

V11.3:

$VersionNumber

Exp[-1000.0]
Precision[%]

We have to explicitly enforce ArbitraryPrecision to work with numbers less than $MinMachineNumber

N[Exp[-1000],10]

Formal definition of Precision[]

A number x with precision p implies that x lies between x(1±((10^(-p))/2))
Precision represents relative uncertainty in representation.

Turns out that a precision of p implies that the first p digits can be trusted.

N[50/26,5]

{50/26(1-0.5*10^-5),50/26 (1+0.5*10^-5)}

SetPrecision[1.2345678,3]

{1.2345678 (1-0.5*10^-3),1.2345678 (1+0.5*10^-3)}

N[ __ ,p] adaptively works towards giving the result with requested precision

Accuracy[]

Similar concept as Precision[] - can be used to enforce ArbitraryPrecision computation
Represents absolute uncertainty. A number x with accuracy a implies that x lies between x±((10^(-a))/2)

Specifying accuracy:

SetAccuracy[1.23487238,10]
Accuracy[%]

1.23487238``10
Accuracy[%]

N[expr, {Infinity,a}] evaluates expr with accuracy a

N[Sqrt[3],{Infinity,10}]
Accuracy[%]

ArbitraryPrecision computation using Accuracy[]:

2-3*(5.2917721067*10^-9)/(.9*10^-13)+3*((5.2917721067*10^-9)/(.9*10^-13))^3-3*(5.2917721067*10^-9)*((.9*10^-13)+(5.2917721067*10^-9))^2*Exp[-2 (.9*10^-13)/(5.2917721067*10^-9)]/(.9*10^-13)^3 (* MachinePrecision computation *)

x=2-3*((Rationalize[5.2917721067*10^-9,0])/(Rationalize[.9*10^-13,0]))+3*((Rationalize[5.2917721067*10^-9,0])/(Rationalize[.9*10^-13,0]))^3-3*Rationalize[5.2917721067*10^-9,0]*(Rationalize[.9*10^-13,0]+Rationalize[5.2917721067*10^-9,0])^2*Exp[-2 (Rationalize[.9*10^-13,0])/(Rationalize[5.2917721067*10^-9,0])]/(Rationalize[.9*10^-13,0])^3
N[x,{Infinity,20}] (* ArbitraryPrecision computation *)

Helpful while dealing with numbers close to zero

Precision or relative uncertainty is undefined for zero

Example 1

poly=1+54 #1+855 #1^2+6300 #1^3+37296 #1^4+221706 #1^5+385782 #1^6+1899783 #1^7-713709 #1^8-8772909 #1^9+1718343 #1^10+17264169 #1^11-3659847 #1^12-20026899 #1^13+3423276 #1^14+13768320 #1^15-1610712 #1^16-5314050 #1^17+881651 #1^18+1545093 #1^19-151263 #1^20-298626 #1^21-24552 #1^22+21978 #1^23+6594 #1^24+792 #1^25+45 #1^26+#1^27&;

r1=Root[poly,1];
Plot[poly[t],{t,r1-.003,r1+.003},ImageSize->Large]

N[poly[r1]]

Mathematica struggles to find the root with requested precision

N[poly[r1],10]
Block[{$MaxExtraPrecision=100000},N[poly[r1],10]]

Works well when accuracy is requested instead

N[poly[r1],{Infinity,10}]

Example 2

N[(Sqrt[3+2 (Sqrt[2])])-1-(Sqrt[2]),30]

(1+(Sqrt[2]))^2//Expand

N[(Sqrt[3+2 (Sqrt[2])])-1-(Sqrt[2]),{Infinity,30}]

Concluding remarks

We learned about:

Three different computation systems used by Mathematica - procedures for setting those up

MachinePrecision is faster than ArbitraryPrecision, however error is not tracked.
With computations involving very small/large numbers we should always re-check with ArbitraryPrecision

Various helper functions for controlling the precision of the results