1

I'm trying to integrate an output variable from a very complex algorithm.

I extensively checked my implementation of the two integration algorithms.

When I use the two methods for my problem, I see a discrepancy: the accuracy is 10-6 and the results are:

Simpson: -1.16844182693551
Romberg: -1.16842500337272

The difference is 1.68 · 10-5, which is greater than the desired accuracy.

Is there any way to know which result is the best?

Cristiano
  • 143
  • 1
    Usually, Romberg is more accurate as the degree of exactness is higher. So without further informations, I would conisder the Romberg-value to be the correct one. – Peter Jun 21 '20 at 09:44
  • @Peter Thank you. Is there anything I can add for a better estimate? – Cristiano Jun 21 '20 at 09:50
  • 1
    Both a description of the problem and details about the application of Simpson and Romberg will be very helpful. Maybe, you just refine your intervals and increase the number of steps in Romberg to get even better results even if this is time-consuming. – Peter Jun 21 '20 at 09:55
  • 1
    What is the accumulated floating point error and/or any other error in your "very complex algorithm"? The higher order columns in the Romberg table will amplify these errors, especially if these are quasi-random and not systematic. – Lutz Lehmann Jun 21 '20 at 09:55
  • @LutzLehmann A good point ! Do you think Simpson is numerically more stable than Romberg ? – Peter Jun 21 '20 at 09:57
  • 1
    @Peter Romberg always stops before reaching the maximum number of steps, in other words, the desired accuracy is always achieved. – Cristiano Jun 21 '20 at 09:58
  • 1
    @LutzLehmann how do I check the accumulated error? I use the algorithm shown here: https://en.wikipedia.org/wiki/Romberg%27s_method#Implementation – Cristiano Jun 21 '20 at 10:01
  • 3
    @Peter : It is a balance, it really depends on the errors involved. Not to forget that Simpson is the 2nd column of the Romberg table, if the midpoint sums are the zeroth column and the trapezoidal values the first column. The lower order methods are more stable for step sizes in their "working range", but require much more evaluations for the same accuracy as higher columns. The higher order columns have a much smaller working range, but get a much better accuracy inside. – Lutz Lehmann Jun 21 '20 at 10:03
  • @Peter It's an astrodynamics problem: first I calculate the position and velocity of an Earth's satellite (it's an algorithm with dozens of trigonometric function), then I calculate the air density at that position (many functions also in that algorithm), finally I calculate the time derivative of the variable that I'm trying to integrate. – Cristiano Jun 21 '20 at 10:05
  • 2
    You can't really check for the accumulated floating point error. The extrapolated error that controls the order is based on the theoretical error from the order of the method. One can estimate where the balance between theoretical error and floating point error will fall, it is a good idea to target this balance. – Lutz Lehmann Jun 21 '20 at 10:06
  • @LutzLehmann the accumulated floating point error should be almost negligible (there are 2 or 3 iterations for the calculation of the position and a few iterations for the air density). For "any other error" I don't know. – Cristiano Jun 21 '20 at 10:15
  • 2
    Much also depends on the function you integrate. An order $p$ formula requires that the function is at least $p$ times differentiable, even more if you consider the error estimation via extrapolation. So in broad strokes, if your model is all smooth functions, then the Romberg value is the more exact one, but its error estimate does not contain the difference from the model to the physical reality. If your model contains actual measurements, then a lower order method should be preferred, but now the estimated error partially also accounts for the distance to physical reality. – Lutz Lehmann Jun 21 '20 at 10:20
  • @LutzLehmann I wouldn't be surprised to know that the air density algorithm generates discontinuities. Neither model (SGP4 and NRLMSISE-00) contains actual measurements. – Cristiano Jun 21 '20 at 10:35

1 Answers1

1

No, given the available information it is not possible to determine which of the two results is the most accurate one. Indeed there is reason to suspect that the error estimate used internally by your routine is faulty or you are outside the range of stepsizes where it is reliable.

However, given a sequence of approximations, say, $A_h, A_{2h}, A_{4h}, \dotsc$ computed using a fixed method and stepsize $h, 2h, 4h, \dotsc,$ we can determine the (effective) order of the method, estimates of the truncation error, as well as the range of stepsizes $h$ for which the error estimates are accurate and the truncation error is larger than the accumulated rounding error. The technique is known as Richardson extrapolation.

These answers to other questions discuss Richardson extrapolation in the context of differentiation, integration and solving ordinary differential equations. The last answer provide the best illustration of what you can do if you compute multiple approximations as it contains actual numerical results.

It is possible to say much more about this, but this is perhaps a subject for another question?

Carl Christian
  • 12,583
  • 1
  • 14
  • 37
  • I tried the Richardson extrapolation method to calculate the derivative of elementary functions and it actually works very well (I didn't know, thank you). Now I try it for the integration. – Cristiano Jun 22 '20 at 08:17