Is the a precision loss due to export/import of numerical values?

Questions
, , ,
1

Hello,

I have the following problem with understanding whether an unwanted
precision loss occurs in my calculations, and I would appreciate some
help. I solve some ODEs using NDSolve, and since the calculations
are very expensive, they are split into several runs. In one run
the calculations are stopped when the available time or memory limits
are approached. The value of the independent variable th, corresponding
to this situation is stored as thmax. This value then becomes the end of the integration interval. The code is:

wprec=50;

ivpsol=NDSolve[
{ivp, WhenEvent[(TimeUsed[] > maxctime || MemoryInUse[] > maxmemory) && th>0, thmax=th;"StopIntegration","LocationMethod"->"StepBegin"]}
,soltable,{th,th0,thmax},WorkingPrecision->wprec,PrecisionGoal->wprec/2,
MaxSteps->Infinity,InterpolationOrder->All
,StepMonitor:>Print[{th,N[Subscript[psix,1][th],20]}]
];

Please note that the working precision is 50, so that one expects
that 50 digits of thmax should be available. Indeed, if I now
export thmax to a file, according to:

Export[file,Join[{thmax},endpsivalues,endphivalues],"List"];

I can see that 50 digits of thmax are in the file (this is
the first of the following numbers stored:

4.2381878739099146976942608527708919906248917152657e6
0.0048230113561383041835207418808668133329086696535244
0.0048249194434038473837547000974380792009777840867437
0.0048330661478961510678031863401831999702176813779798
etc.

I also print the value into another file, according to:

Write[strm,OutputForm["# thmax = "],CForm[thmax]];

with the result:

# thmax = 4.2381878739099146976942608527708919906248917152657e6

However, in a next run of my calculations, the file “file” is imported
to provide starting values for further integration, and the former variable
thmax becomes the initial value th0 of the independent variable,
according to the code:

initvalues=SetPrecision[Import[file,"List"],Infinity];
th0=First[initvalues];

When I now print the th0, according to the code:

Write[strm,OutputForm["# th0 = "],CForm[th0]];

I get the printout:

# th0 = 4.238187873909915e6

As it seems, only 16 digits are printed, which seems to indicate some
precision loss, despite the fact that SetPrecision[…,Infinity] was used
on input.

Is this really a precision loss, and how can one avoid it?

Lesław

2

I am no expert in numerics, and I am just guessing. When you SetPrecision[…,Infinity], you are converting the floating point number into an exact rational number. When you insert the rational number back into NDSolve, it gets converted back to floating point, according to some criterion, possibly MachinePrecision. For example, starting with all exact numbers you get MachinePrecision domain:

sol = NDSolveValue[{x'[t] == x[t], x[0] == 1}, x, {t, 0, 1}]
sol["Domain"]
Precision[%]
3

In the code I send I did not go up to this point. The printout of th0 is independent of what happens next with this value.

But, I can add that in NDSolve I use WorkingPrecision=50, so that in my opinion all input to NDSolve, including th0, should be converted to this precision (from Infinite precision, intentionally). If th0 was of lower precision, I would expect NDSolve to complain about insufficient precision of input, but it does not complain.

Lesław

4

Even if you use a high WorkingPrecision, in the calculation the actual output precision drops.

5

If I copy and paste the numbers you posted as shown below, the precision of the numbers is slightly less than 50, from 49.6 to 49.7, due to conversion from the ASCII representation:

{4.2381878739099146976942608527708919906248917152657*^6, 
0.0048230113561383041835207418808668133329086696535244,
 0.0048249194434038473837547000974380792009777840867437, 
0.0048330661478961510678031863401831999702176813779798}

This has nothing to do with the main issue, but I felt it important to realize what numbers I am working with. If I export and import them, there is a further slight loss in precision, again due to the round trip conversion to and from ASCII. Again, this is not particularly important, since the precision is 47.6 to 48.7, a loss of 1-2 digits. It might be important to realize that in the conversion to ASCII, trailing zeros are omitted, and this translates to an extra precision loss equal to the number of zeros dropped. Whether this should be a significant concern depends on how the numbers are used in further calculations.

If I promote the precision of the numbers to 50 by adding a tick-50 at the ends, then the export/import loss of precision is the same as the copy-paste one, about 0.3–0.4 digits.

Export[FileNameJoin[{$TemporaryDirectory, "foo"}],
 {4.2381878739099146976942608527708919906248917152657`50*^6, 
  0.0048230113561383041835207418808668133329086696535244`50,
  0.0048249194434038473837547000974380792009777840867437`50, 
  0.0048330661478961510678031863401831999702176813779798`50}, "List"]
Import[FileNameJoin[{$TemporaryDirectory, "foo"}], "List"] // FullForm

In short, the data has only a minor loss of precision, nothing near the 34-digit loss you see in CForm[th0]. Let’s look at CForm[]. First, th0 is a rational number (infinite precision). CForm[th0] is the CForm of a rational number (check with FullForm below). What you’re seeing is the effect of typesetting the CForm of a rational number by the Front End. Try CForm[1/3] for a simpler example. For whatever reason (I did not find it mentioned in the documentation), what you get is a C double approximation of the rational number (update: SE Q&A on this topic). The upshot is that CForm[] is misleading you as to the precision of th0.

initvalues = SetPrecision[Import[FileNameJoin[{$TemporaryDirectory, "foo"}], "List"], Infinity];
th0 = First[initvalues]
(*Out[]=   1476792935618238128801055858742073436930432589667/348449143727040986586495598010130648530944  *)

CForm[th0] (* defaults to machine-precision form *)
(*Out[]=  4.238187873909915e6  *)

CForm[N[th0, 50]]
(*Out[]=  4.2381878739099146976942608527708919906248917152657e6  *)

CForm[th0] // FullForm
(*Out=  CForm[Rational[
  1476792935618238128801055858742073436930432589667,
  348449143727040986586495598010130648530944]]  *)

You might want to consider exporting the WL form of the numbers you wish to save. These are imported in the same form as they were exported. Of course, the files cannot be catenated, but I suppose that is unimportant.

Export[FileNameJoin[{$TemporaryDirectory, "foo"}],
 {4.2381878739099146976942608527708919906248917152657`50*^6, 
  0.0048230113561383041835207418808668133329086696535244`50,
  0.0048249194434038473837547000974380792009777840867437`50, 
  0.0048330661478961510678031863401831999702176813779798`50}, "WL"]
Import[FileNameJoin[{$TemporaryDirectory, "foo"}], "WL"] // FullForm

Alternatively, you might wish to store the extra guard digits. They will import with a higher precision than 50 (unless there are trailing zeros), but you could use SetPrecision[] to set the precision back to 50.

x = N[Pi, 50];
Export[FileNameJoin[{$TemporaryDirectory, "foo"}],
 {First@StringSplit[ToString[x, InputForm], "`"]}, "List"]
Import[FileNameJoin[{$TemporaryDirectory, "foo"}], "List"] // FullForm