# -*- coding: utf-8 -*-
# mathtoolspy
# -----------
# A fast, efficient Python library for mathematically operations, like
# integration, solver, distributions and other useful functions.
#
# Author: sonntagsgesicht, based on a fork of Deutsche Postbank [pbrisk]
# Version: 0.3, copyright Wednesday, 18 September 2019
# Website: https://github.com/sonntagsgesicht/mathtoolspy
# License: Apache License 2.0 (see LICENSE file)
[docs]class GaussLobattoStep:
def __init__(self, function, alpha, beta, Is, max_number_of_iterations):
self._function = function
self._alpha = alpha
self._beta = beta
self._is = Is
self._max_number_of_iterations = max_number_of_iterations
self._used_iter = 0
[docs] def adaptive_step(self, a, b, fa, fb):
h = (b - a) / 2.0
m = (a + b) / 2.0
mll = m - self._alpha * h
ml = m - self._beta * h
mr = m + self._beta * h
mrr = m + self._alpha * h
fmll = self._function(mll)
fml = self._function(ml)
fm = self._function(m)
fmr = self._function(mr)
fmrr = self._function(mrr)
i2 = (h / 6.0) * (fa + fb + 5.0 * (fml + fmr))
i1 = (h / 1470.0) * (77.0 * (fa + fb) + 432.0 * (fmll + fmrr) + 625.0 * (fml + fmr) + 672.0 * fm)
if self._is + (i1 - i2) == self._is or mll <= a or b <= mrr:
if m <= a or b <= m:
raise Exception("Interval contains no more machine number.")
return i1
self._used_iter += 1
if self._used_iter > self._max_number_of_iterations:
raise Exception("Maximal number of iterations " + str(self._max_number_of_iterations) + "'exceeded.")
ret1 = self.adaptive_step(a, mll, fa, fmll)
ret2 = self.adaptive_step(mll, ml, fmll, fml)
ret3 = self.adaptive_step(ml, m, fml, fm)
ret4 = self.adaptive_step(m, mr, fm, fmr)
ret5 = self.adaptive_step(mr, mrr, fmr, fmrr)
ret6 = self.adaptive_step(mrr, b, fmrr, fb)
return ret1 + ret2 + ret3 + ret4 + ret5 + ret6
[docs]class GaussLobattoIntegrator:
def __init__(self, max_number_of_iterations=255, abs_tolerance=1.0e-10):
self.max_number_of_iterations = max_number_of_iterations
self.abs_tolerance = abs_tolerance
[docs] def integrate(self, function, lower_bound, upper_bound):
if lower_bound == upper_bound:
return 0
if lower_bound > upper_bound:
return -self.integrate(function, upper_bound, lower_bound)
# three constants for the Gauss-Lobatto approach
x1 = 0.942882415695480
x2 = 0.641853342345781
x3 = 0.236383199662150
alpha = 0.816496580927726 # sqrt(2.0/3.0)
beta = 0.447213595499958 # 1.0 / sqrt(5.0)
m = (lower_bound + upper_bound) / 2.0
h = (upper_bound - lower_bound) / 2.0
fa = function(lower_bound)
fb = function(upper_bound)
f2 = function(m - x1 * h)
f3 = function(m - alpha * h)
f4 = function(m - x2 * h)
f5 = function(m - beta * h)
f6 = function(m - x3 * h)
f7 = function(m)
f8 = function(m + x3 * h)
f9 = function(m + beta * h)
f10 = function(m + x2 * h)
f11 = function(m + alpha * h)
f12 = function(m + x1 * h)
w1 = 0.158271919734802
w2 = 0.0942738402188500
w3 = 0.155071987336585
w4 = 0.188821573960182
w5 = 0.199773405226859
w6 = 0.224926465333340
w7 = 0.242611071901408
i2 = (h / 6.0) * (fa + fb + 5.0 * (f5 + f9))
i1 = (h / 1470.0) * (77.0 * (fa + fb) + 432.0 * (f3 + f11) + 625.0 * (f5 + f9) + 672.0 * f7)
Is = h * (w1 * (fa + fb)) + w2 * (f2 + f12) + w3 * (f3 + f11) + w4 * (f4 + f10) + w5 * (f5 + f9) + w6 * (
f6 + f8) + w7 * f7
sign = 1 if Is >= 0 else -1
R = abs(i1 - Is) / abs(i2 - Is)
tol = self.abs_tolerance
tol = tol / R if R > 0 and R < 1 else tol
Is = sign * abs(Is) * tol / 1E-12
Is = upper_bound - lower_bound if Is == 0 else Is
gls = GaussLobattoStep(function, alpha, beta, Is, self.max_number_of_iterations)
ret = gls.adaptive_step(lower_bound, upper_bound, fa, fb)
return ret
def __call__(self, fct, lower_bound, upper_bound):
return self.integrate(fct, lower_bound, upper_bound)