Please, consider the following MWE representing my code:

analytic = y /. DSolve[{y'[x] == y[x], y[0] == 1}, y, x][[1]];
numeric = y /. NDSolve[{y'[x] == y[x], y[0] == 1}, y, {x, 0, 1}, WorkingPrecision -> precision][[1]];
numeric >> "test.m";
interpolated = Interpolation[Table[{x, numeric[x]}, {x, Range[0, 1, 1/100]}]];
loaded = << "test.m";

Let's evaluate the expressions

{numeric[%], interpolated[%], loaded[%]};
Map[Precision, %]
%% - analytic[%%%]

For precision = 20 we get

{19.4718, 19.4718, MachinePrecision}
{-1.82359038525*10^-8, -1.82359038525*10^-8, -1.82359*10^-8}

For precision = 40 we get

{39.7727, 39.7727, MachinePrecision}
{-3.4833827190165031850763300823*10^-11, -3.4833827190165031850763300823*10^-11, -3.48339*10^-11}

First, while the values in the file test.m are saved with the given precision, when loaded back into the notebook, only MachinePrecision values are used. How should I (save and) load the results so that the precision is kept?

I do understand that the solution is going to have lower precision than the one given to WorkingPrecision. Nevertheless, the WorkingPrecision have doubled but the error of the solution is only slightly smaller and well below the MachinePrecision.

Moreover, when precision is set to the default value MachinePrecision the result is

{MachinePrecision, MachinePrecision, MachinePrecision}
{-1.24619*10^-9, -1.24619*10^-9, -1.24619*10^-9}

While in my case MachinePrecision = 15.9546 the result is better than in case of precision = 20. Why is that?

I would like to have some kind of precision control such that doubling a parameter at the start of the computation would roughly translate in double the amount of correct digits at the end of the computation. How can I achieve this? All user given numbers are of infinite precision.

(I run the code inside ParallelSubmit using wolframscript.)


The Stack Exchange question pointed by @Domen solves the issue with loading the solution. Nevertheless, the other more important issues still remain.

Of course, I did set the particular precision for the corresponding version of numeric, I just didn't spell out the full code for brevity.

I am not surprised to lose precision when subtracting two close numbers. I am surprised by how little is gained by starting with double the precision a even more so why the lower MachinePrecision gives better result than precision = 20.

  • 4
    $\begingroup$ Does this answer your question? Precision of interpolation function after exporting and importing $\endgroup$
    – Domen
    Commented Apr 23 at 20:48
  • $\begingroup$ @Domen Thank you for pointing out this question. It does help with part of the question, the main part of the question is, however, about why the precision (outside put and get) behaves as discussed and how to make to work as intended. $\endgroup$
    – fales
    Commented Apr 24 at 17:45

2 Answers 2


Ignoring the issue solved in the linked question (proposed duplicate), the extra error is due to interpolation error. The default interpolation order in a solution to a $k$th-order ODE returned by NDSolve is $2k+1$, or $3$ in the OP's example. The actual order of the solution is usually much higher and varies in the default multistep LSODA method used by NDSolve. To get full accuracy, use InterpolationOrder -> All.

For precision = 40, the error is about 10^-20, which accords with the default PrecisionGoal of WorkingPrecision/2:

precision = 40;
analytic = y /. DSolve[{y'[x] == y[x], y[0] == 1}, y, x][[1]];
numeric = 
  y /. NDSolve[{y'[x] == y[x], y[0] == 1}, y, {x, 0, 1}, 
     WorkingPrecision -> precision, InterpolationOrder -> All][[1]];
interpolated = 
  Interpolation[Table[{x, numeric[x]}, {x, Range[0, 1, 1/100]}]];
{numeric[%], interpolated[%]} - analytic[%]

(*  {5.3733064813649050869*10^-20, 5.3733064813649050869*10^-20}  *)

Another way to look at it is that the error at the steps computed by NDSolve are all on the order of 5 * 10^-20 or less, whether InterpolationOrder -> All is included or omitted. Both ways yield the same result below:

xgrid = numeric["Grid"] // Flatten;
With[{exact = analytic[xgrid]},
 (numeric[xgrid] - exact)/(10^-precision + Abs@exact) // Abs // Max]

(*  5.9160946879896551988*10^-20  *)

The error using interpolated[xgrid] is on the order of 5±1 * 10^-10, using the two methods. You still get a large error because the steps in xgrid do not coincide with the interpolation grid for interpolated[], and you get extra error from the interpolation when interpolated[xgrid] is evaluated.


(* "14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023)" *)


analytic = y /. DSolve[{y'[x] == y[x], y[0] == 1}, y, x][[1]];

If you define numeric using Set, it will be permanently defined using machine precision. Use SetDelayed with the precision as a parameter.

numeric[prec_?NumericQ] := y /. NDSolve[
     {y'[x] == y[x], y[0] == 1}, y, {x, 0, 1},
     WorkingPrecision -> prec][[1]];

(analytic[x] - numeric[20][x]) /. x -> 11/100

(* 1.82359038525*10^-8 *)

The loss in precision is

20 - Precision[%]

(* 8.31507 *)

(analytic[x] - numeric[40][x]) /. x -> 11/100

(* 3.4833827190165031850763300823*10^-11 *)

The loss in precision is

40 - Precision[%]

(* 10.7331 *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.