Day 98: Romberg integration
Romberg’s method to find a definite integral combines two formulas, extended trapezoidal rule and Richardson extrapolation, to get a good approximation in relatively low number of steps.
Let’s explain how it all works without causing a serious headache.
The simplest way to find a definite integral of function f on interval (a, b) is to use a trapezoidal rule.
It’s the left-most chart and you can see that the formula is merely an area of trapezoid.
You can also notice that there’s a certain error which can be improved using extended trapezoidal rule. Split interval (a, b) into two halves, apply trapezoidal rule onto each half and sum them up.
It’s the chart in the middle. The formula is merely an area of two trapezoids.
The third plot contains another halving and even better estimate. We can go on and on until we are satisfied with the estimate error.
The problem is, the convergence in this way would be very slow. And here comes Richardson extrapolation to the rescue.
In Richardson extrapolation we take two bad estimates and combine them into a good one. The best explanation I’ve ever seen is definitely described in this Dilbert strip, check it out.
However, unlike Dilbert we are lucky since we know error bound of each estimate.
In other words, if we are sure that one estimate is better that the other [e.g. due to smaller interval in extended trapezoidal rule], we also know in which direction should we expect the exact solution.
Finally, Romberg’s method builds a table of estimates with surprisingly accurate results.
The first column contains results of extended trapezoidal rule, each row using twice as many splits as the row before. All the other values are results of Richardson extrapolation.
The diagonal contains the final estimates and each diagonal element gives about two more digits of precision. Hence the element at 5th row and 5th column gives an estimate of about 10 digits of accuracy.
https://github.com/coells/100days
https://notebooks.azure.com/coells/libraries/100days
algorithm
def integrate(fn, a, b, steps=5, debug=False, exact=None):
table = np.zeros((steps, steps), dtype=np.float64)
pow_4 = 4 ** np.arange(steps, dtype=np.float64) - 1 # trapezoidal rule
h = (b - a)
table[0, 0] = h * (fn(a) + fn(b)) / 2 for j in range(1, steps):
h /= 2 # extended trapezoidal rule
table[j, 0] = table[j - 1, 0] / 2
table[j, 0] += h * np.sum(
fn(a + i * h) for i in range(1, 2 ** j + 1, 2)
) # richardson extrapolation
for k in range(1, j + 1):
table[j, k] = table[j, k - 1] + \
(table[j, k - 1] - table[j - 1, k - 1]) / pow_4[k] # debug
if debug:
print(table, file=sys.stderr)
if exact is not None:
errors = [
'%.2e' % i
for i in np.abs(table.diagonal() - exact)
]
print('abs. error:', errors, file=sys.stderr) return table[-1, -1]
integral[0, 1] of e^(-x²)
integrate(lambda x: np.exp(-x * x), 0, 1,
debug=True, exact=0.746824132812427)0.74682401848228175[[ 0.6839397206 0. 0. 0. ]
[ 0.7313702518 0.7471804289 0. 0. ]
[ 0.7429840978 0.7468553798 0.7468337098 0. ]
[ 0.7458656148 0.7468261205 0.7468241699 0.7468240185]]abs. error: ['6.29e-02', '3.56e-04', '9.58e-06', '1.14e-07']
ln(2)
integrate(1..__truediv__, 1, 2, debug=True, exact=np.log(2))0.6931474776448322[[ 0.75 0. 0. 0. ]
[ 0.7083333333 0.6944444444 0. 0. ]
[ 0.6970238095 0.6932539683 0.6931746032 0. ]
[ 0.6941218504 0.6931545307 0.6931479015 0.6931474776]]abs. error: ['5.69e-02', '1.30e-03', '2.74e-05', '2.97e-07']
integral[0, 1] of x³
integrate(lambda x: x**3, 0, 1, debug=True, exact=.25)0.25[[ 0.5 0. 0. 0. ]
[ 0.3125 0.25 0. 0. ]
[ 0.265625 0.25 0.25 0. ]
[ 0.25390625 0.25 0.25 0.25 ]]abs. error: ['2.50e-01', '0.00e+00', '0.00e+00', '0.00e+00']
logarithmus naturalis
def ln(x):
if x <= 0:
raise ValueError()
m, e = np.frexp(x)
return integrate(1..__truediv__, 1, m) + e * 0.6931471805599453
ln(np.e)0.99999999844595155ln(2)0.69314688347505837ln(np.pi)1.1447298858266319
normal distribution
def norm_pdf(x, mean, sd):
x0 = (x - mean) ** 2
v2 = 2 * sd ** 2
return np.exp(-x0 / v2) / np.sqrt(np.pi * v2)def norm_cdf(x, mean, sd):
return integrate(lambda x: norm_pdf(x, mean, sd), mean, x) + .5
norm_cdf(1.96, mean=0, sd=1)0.975000070023468norm_cdf(1.2, mean=1, sd=.2) - norm_cdf(0.8, mean=1, sd=.2)0.68268949398919654
