Integration Romberg Method with variable step

The Romberg Methods is one of the most popular algorithm for numeric integration of analytical functions. This method has several very attractive features:

• It is easy to translates it by any programming language.
• It is a very robust method.
• It is extremely fast. On 32-bit, it can reach the machine precision in only 3-4 extrapolation.
• It is also very efficients: usually it requires less functions computing (nodes) then other methods.
• It do not need different complicated formulas: The most part of computing is done by the simple arithmetic addition. Only (few) extrapolations require also the division.

In the example below we compare integration result with different method

The first 5 curve are obtain by Newton-Cotes formulas with different degree (Degree= 1 is the Trapezoidal rule, Degree=2 is the Cavalieri-Simpson formula, Degree=4 is the Bode formula)

As we can see, Romberg integration reach the precision of about 5E-14 in 512 nodes, and it is comparable with a 6° degree integration polynomial. It is a remarkable goal but can it do better? The answer is yes, at least for this type of function.

Splitting the integration interval between two sub-intervals [1 , 3] and [3 ,. 10], we get:

This simple trik is the basic of step variable integration. The reason of the saving computational effort must be searched in the fact that the truncation error is very disomogeneous along all the integration interval.

Derivative of 2° order

If the step h is fixed the error estimation can show by the 2° derivative of integration function f(x). Using the same function samples of the integration, the second derivative is approximated by the following interpolation formula:

The following graph shows the curve obtained from 64 nodes of given function.

As we ca see that for that in near the initial value x=1 the f'''(x) , and also the truncation error, is very higher than the one at the final point x=10. So the integral over the rigth-subinterval [3 , 10] can be calculated with a step higher than the left one. This simple trick save several computation values of f(x).

Splitting Rule - The Romberg method with variable step uses the 2° derivative formula to extimate the error truncation along the interval. It does not involve any more computation of f(x) because it use only the same samples of integration.

For the same reason the interval is always divided in two half part (splitting rule).

|-----|-----|-----|-----|-----|-----|-----|-----|

a    a+h   a+2h  a+3h  a+4h  a+5h  a+6h  a+7h   b

For example, if the interval [a , b] is divided in 8 part, the we split it in two sub-interval at the middle point a+4h. So we can campute two separate integral into the interval [a, a+4h] [a+4h , b]

Note that for each interval we can use the previsous sample value of f(x).

Splitting factor - To decide when it splits an integration interval, the algorithm has to make same calculation. One of the simpler is the following factor:

Where f'' are the 2° derivative values obtained from the above given formula.

The process are recursively applyed to the new intervals: each split gives two new intervals with half range, and so on. The number of split depends on the regular of function. Usually a limit of max number of split is set (32, 64)

For each sub-intervals a Romberg integration is performed, abtaining separate integral values and error values. If an error of a sub-interva falls under a the specific Emax, then the computation is stopped only for this sub-interval. When all the sub-intervals are stopped, the process ends. The final integral is the summation of all sub-integrals.

A routine that realized this algorithm is in the workbooks

Applying this method to the initial test integral we have: 0.900000000000078 (Error = -7.7716E-14 , 225 nodes)

having splitted the original [1 , 10] interval in the following sub-intervals

The standard Romberg integration along the same interval has given about 5E-14 precision with 512 nodes while here we have get the simila precision with only 225 nodes; so the efficency has increased of more then 60%.

The efficency increases much more for larger integration interval, (see: Results of variable step intgration):

HOME