]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Mention Kahan summation and discuss 'double' vs 'long double'. 13200/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Fri, 7 Jan 2022 19:52:12 +0000 (12:52 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Sat, 8 Jan 2022 20:04:43 +0000 (13:04 -0700)
examples/step-10/doc/results.dox
examples/step-10/step-10.cc

index 628bfc69d20a6d74da38ae51010d75f082a6c0cb..8b6e7a0471752479216853d5caeccb0345e10cb4 100644 (file)
@@ -190,7 +190,7 @@ cells      eval.pi            error
   320 3.1415926535897576 3.5527e-14 7.94
  1280 3.1415926535897896 3.5527e-15 3.32
  5120 3.1415926535897940 8.8818e-16 2.00
- @endcode
+@endcode
 
 @note Once the error reaches a level on the
   order of $10^{-13}$ to $10^{-15}$, it is essentially dominated by
@@ -198,6 +198,9 @@ cells      eval.pi            error
   in internal computations. Since these things change, the precise values
   and errors change from release to release at these round-off levels,
   though the overall order of errors should of course remain the same.
+  See also the comment below in the section on
+  <a href="#extensions">Possibilities for extensions</a> about how to compute
+  these results more accurately.
 
 One of the immediate observations from the output above is that in all cases
 the values converge quickly to the true value of
@@ -223,3 +226,97 @@ convergence. (This effect is also discussed in detail in the following
 publication: A. Bonito, A. Demlow, and J. Owen: "A priori error
 estimates for finite element approximations to eigenvalues and
 eigenfunctions of the Laplace-Beltrami operator", submitted, 2018.)
+
+
+
+<h3>Possibilities for extensions</h3>
+<a name="extensions"></a>
+
+As the table of numbers copied from the output of the program shows above,
+it is not very difficult to compute the value of $\pi$ to 13 or 15 digits. But,
+the output also shows that once we approach the level of accuracy with which
+`double` precision numbers store information (namely, with roughly 16 digits
+of accuracy), we no longer see the expected convergence order and the error
+no longer decreases with mesh refinement as anticipated. This is because both
+within this code and within the many computations that happen within deal.II
+itself, each operation incurs an error on the order of $10^{-16}$; adding
+such errors many times over then results in an error that may be on the
+order of $10^{-14}$, which will dominate the discretization error after
+a number of refinement steps and consequently destroy the convergence rate.
+
+The question is whether one can do anything about this. One thought is to
+use a higher-precision data type. For example, one could think of declaring
+both the `area` and `perimeter` variables in `compute_pi_by_area()` and
+`compute_pi_by_perimeter()` with data type `long double`. `long double`
+is a data type that is not well specified in the C++ standard but at least
+on Intel processors has around 19, instead of around 16, digits of accuracy.
+If we were to do that, we would get results that differ from the ones shown
+above. However, maybe counter-intuitively, they are not uniformly better.
+For example, when computing $\pi$ by the area, at the time of writing
+these sentences we get these values with `double` precision for degree 4:
+@code
+    5 3.1415871927401144 5.4608e-06    -
+   20 3.1415926314742491 2.2116e-08 7.95
+   80 3.1415926535026268 8.7166e-11 7.99
+  320 3.1415926535894005 3.9257e-13 7.79
+ 1280 3.1415926535899774 1.8430e-13 1.09
+ 5120 3.1415926535897669 2.6201e-14 2.81
+@endcode
+On the other hand, the results are as follows when using `long double`:
+@code
+ cells      eval.pi            error
+    5 3.1415871927401136 5.4608e-06    -
+   20 3.1415926314742446 2.2116e-08 7.95
+   80 3.1415926535026215 8.7172e-11 7.99
+  320 3.1415926535894516 3.4157e-13 8.00
+ 1280 3.1415926535897918 1.5339e-15 7.80
+ 5120 3.1415926535897927 5.2649e-16 1.54
+@endcode
+Indeed, here we get results that are approximately 50 times as accurate.
+On the other hand, when computing $\pi$ by the perimeter, we get this with
+`double` precision:
+@code
+    5 3.1415921029432572 5.5065e-07    -
+   20 3.1415926513737582 2.2160e-09 7.96
+   80 3.1415926535810699 8.7232e-12 7.99
+  320 3.1415926535897576 3.5527e-14 7.94
+ 1280 3.1415926535897896 3.5527e-15 3.32
+ 5120 3.1415926535897940 8.8818e-16 2.00
+@endcode
+Whereas we get the following with `long double`:
+@code
+    5 3.1415921029432572 5.5065e-07     -
+   20 3.1415926513737595 2.2160e-09  7.96
+   80 3.1415926535810703 8.7230e-12  7.99
+  320 3.1415926535897576 3.5705e-14  7.93
+ 1280 3.1415926535897918 1.3785e-15  4.70
+ 5120 3.1415926535897944 1.3798e-15 -0.00
+@endcode
+Here, using `double` precision is more accurate by about a factor of
+two. (Of course, in all cases, we have computed $\pi$ with more
+accuracy than any engineer would ever want to know.)
+
+What explains this unpredictability? In general, round-off errors can
+be thought of as random, and add up in ways that are not worth thinking
+too much about; we should therefore always treat any accuracy beyond, say,
+thirteen digits as suspect. Thus, it is probably not worth spending
+too much time on wondering why we get different winners and losers in the
+data type exchange from `double` and `long double`. The accuracy of the
+results is also largely not determined by the precision of the data type
+in which we accumulate each cell's (or face's) contributions, but the
+accuracy of what deal.II gives us via FEValues::JxW() and FEFaceValues::JxW(),
+which always uses `double` precision and which we cannot directly affect.
+
+But there are cases where one can do something about the precision, and it
+is worth at least mentioning the name of the most well-known algorithm in
+this area. Specifically, what we are doing when we add contributions into
+the `area` and `perimeter` values is that we are adding together *positive*
+numbers as we do here. In general, the round-off errors associated with each
+of these numbers is random, and if we add up contributions of substantially
+different sizes, then we will likely be dominated by the error in the largest
+contributions. One can avoid this by adding up numbers sorted by their
+size, and this may then result in marginally more accurate end results.
+The algorithm that implements this is typically called
+<a href="https://en.wikipedia.org/wiki/Kahan_summation_algorithm">Kahan's summation algorithm</a>.
+While one could play with it in the current context, it is likely not going
+to improve the accuracy in ways that will truly matter.
index 88bb451ee3c3144f96f54b55d49ca83410b48cdf..237207efd906e0467eaf152a85b5c23c561f334e 100644 (file)
@@ -349,9 +349,7 @@ namespace Step10
                   }
             // Then store the evaluated values in the table...
             table.add_value("eval.pi", static_cast<double>(perimeter / 2.0));
-            table.add_value("error",
-                            static_cast<double>(
-                              std::fabs(perimeter / 2.0 - numbers::PI)));
+            table.add_value("error", std::fabs(perimeter / 2.0 - numbers::PI));
           }
 
         // ...and end this function as we did in the previous one:

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.