.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_numerical_methods/general_methods/plot_estimate_integral_Gauss_Kronrod.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_numerical_methods_general_methods_plot_estimate_integral_Gauss_Kronrod.py: Integrate a function with Gauss-Kronrod algorithm ================================================= .. GENERATED FROM PYTHON SOURCE LINES 7-11 .. code-block:: Python import openturns as ot import numpy as np import openturns.viewer as otv .. GENERATED FROM PYTHON SOURCE LINES 12-16 References : - [davis1975]_ section "2.7.1.1 "The Kronrod scheme", p. 82 and section 6.2.0 "An iterative nonadaptive scheme based on Kronrod formulas" p.321. - [dahlquist2008]_ p.575. .. GENERATED FROM PYTHON SOURCE LINES 18-48 Introduction ============ In this example, we present the :class:`~openturns.GaussKronrod` algorithm for one dimensional integration. That is, the algorithm can approximate the integral: .. math:: \int_a^b f(x) dx where :math:`f:[a,b] \rightarrow \mathbb{R}^p` is a function, :math:`[a,b] \subset \mathbb{R}` with :math:`a\leq b` is a one dimensional interval, :math:`p` is the dimension of the output. Notice that the dimension of the input must be equal to 1, but the number of outputs can be greater than 1. Suppose that we have estimated the integral with Gaussian quadrature and :math:`m` quadrature nodes. If we want to improve the accuracy and use more nodes, the issue is that the new nodes do not correspond to the old ones: therefore, we cannot reuse the function evaluations. The Gauss-Kronrod algorithm improves the situation by using two different methods: - a Gaussian quadrature rule with :math:`m` nodes, - a Kronrod extension with :math:`2m+1` nodes. The rule :math:`(G_m,K_{2m+1})` is called a Gauss-Kronrod pair. In the Kronrod extension, the first :math:`m` nodes are equal to the nodes in Gaussian quadrature. The Gaussian quadrature rule with :math:`m` nodes is exact for polynomials of degree :math:`2m-1`. The Kronrod extension with :math:`2m+1` nodes is designed to be exact for polynomials of degree :math:`3m+1`. The choice of the weight function :math:`w(x)` determines the nodes. We consider the weight :math:`w(x)=1` and the interval :math:`[a,b]=[-1,1]` (it is straightforward to generalize this for an arbitrary interval :math:`[a,b]`). In this case, the new :math:`m+1` nodes of the Kronrod extension interlaces with the Gaussian nodes. The weights are guaranteed to be positive (an essential property for numerical stability). .. GENERATED FROM PYTHON SOURCE LINES 50-56 Example ======= The following example is from [davis1975]_ p.325: .. math:: \int_0^1 f(x) dx = \int_0^1 \frac{2}{2 + \sin(10 \pi x)} dx = \frac{2}{\sqrt{3}} = 1.154700538379251529. .. GENERATED FROM PYTHON SOURCE LINES 58-59 We first define the function as a :class:`~openturns.SymbolicFunction`. .. GENERATED FROM PYTHON SOURCE LINES 59-64 .. code-block:: Python integrand = ot.SymbolicFunction(["x"], ["2 / (2 + sin(10 * pi_ * x))"]) integrand.setOutputDescription([r"$\frac{2}{2 + sin(10 \pi x)}$"]) graph = integrand.draw(0.0, 1.0, 200) _ = otv.View(graph) .. image-sg:: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_Gauss_Kronrod_001.png :alt: $\frac{2}{2 + sin(10 \pi x)}$ as a function of x :srcset: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_Gauss_Kronrod_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 65-69 We see that regular spikes and valleys will make this function difficult to integrate, because of the large curvatures implied at these points. We will later count the number of function evaluations. But a small amount of function evaluations has already been used for the plot and this is why we must take it into account. .. GENERATED FROM PYTHON SOURCE LINES 69-72 .. code-block:: Python before_evaluation_number = integrand.getEvaluationCallsNumber() before_evaluation_number .. rst-class:: sphx-glr-script-out .. code-block:: none 200 .. GENERATED FROM PYTHON SOURCE LINES 73-77 Basic use ========= We first choose the Gauss-Kronrod rule. Six quadratures are available: we select the "G11K23" rule. It uses 11 nodes from a Gauss rule and 23 nodes from a Kronrod rule, re-using the nodes from the Gauss rule. .. GENERATED FROM PYTHON SOURCE LINES 77-79 .. code-block:: Python quadrature_rule = ot.GaussKronrodRule(ot.GaussKronrodRule.G11K23) .. GENERATED FROM PYTHON SOURCE LINES 80-81 We set the maximum number of sub-intervals and the maximum absolute error. .. GENERATED FROM PYTHON SOURCE LINES 81-88 .. code-block:: Python maximumSubIntervals = 100 maximumError = 1.0e-8 algo = ot.GaussKronrod(maximumSubIntervals, maximumError, quadrature_rule) interval = ot.Interval(0.0, 1.0) computed = algo.integrate(integrand, interval) computed[0] .. rst-class:: sphx-glr-script-out .. code-block:: none 1.1547005383792528 .. GENERATED FROM PYTHON SOURCE LINES 89-93 Notice that the algorithm can integrate a function which has several outputs (but the number of inputs is restricted to 1). This is why we use the index `[0]` of `computed`, since :meth:`~openturns.GaussKronrod.integrate` returns a :class:`~openturns.Point`. In order to check this computation, we compute the log-relative error in base 10. In most cases (except when the exponent of the two numbers are different), this represents the number of correct digits in base 10. .. GENERATED FROM PYTHON SOURCE LINES 93-97 .. code-block:: Python expected = 1.154700538379251529 LRE_10 = -np.log10(abs(computed[0] - expected) / abs(expected)) LRE_10 .. rst-class:: sphx-glr-script-out .. code-block:: none np.float64(14.937877892447528) .. GENERATED FROM PYTHON SOURCE LINES 98-101 The method computes more than 14 digits correctly. Given that 17 digits is the best we can, this is an astonishing performance. We then compute the number of function evaluations. .. GENERATED FROM PYTHON SOURCE LINES 101-105 .. code-block:: Python after_evaluation_number = integrand.getEvaluationCallsNumber() number_of_calls = after_evaluation_number - before_evaluation_number number_of_calls .. rst-class:: sphx-glr-script-out .. code-block:: none 506 .. GENERATED FROM PYTHON SOURCE LINES 106-113 Advanced use ============ The Gauss-Kronrod algorithm strives to produce an approximated integral which actual error is less than the tolerance. The algorithm estimates the error, which may be used to guess the accuracy in the situation where the exact value is unknown (this is the general use case, of course). In order to get the error estimated by the algorithm, we use the third input argument of the :meth:`~openturns.GaussKronrod.integrate` method. .. GENERATED FROM PYTHON SOURCE LINES 113-117 .. code-block:: Python error = ot.Point() computed = algo.integrate(integrand, interval, error) computed[0] .. rst-class:: sphx-glr-script-out .. code-block:: none 1.1547005383792528 .. GENERATED FROM PYTHON SOURCE LINES 118-119 The variable `error` now contains the error estimate from the algorithm. .. GENERATED FROM PYTHON SOURCE LINES 119-121 .. code-block:: Python error[0] .. rst-class:: sphx-glr-script-out .. code-block:: none 1.6447801423976603e-09 .. GENERATED FROM PYTHON SOURCE LINES 122-123 We see that the error estimate is a little lower than the tolerance, which indicates that the integral should be correctly approximated. .. GENERATED FROM PYTHON SOURCE LINES 125-128 In the next advanced example, let us now use the powerful feature of the algorithm that is, the ability to get the intermediate sub-intervals used by the algorithm. .. GENERATED FROM PYTHON SOURCE LINES 128-138 .. code-block:: Python error = ot.Point() lowerBound = 0.0 upperBound = 1.0 ai = ot.Point() bi = ot.Point() ei = ot.Point() fi = ot.Sample() computed = algo.integrate(integrand, lowerBound, upperBound, error, ai, bi, fi, ei) computed[0] .. rst-class:: sphx-glr-script-out .. code-block:: none 1.1547005383792528 .. GENERATED FROM PYTHON SOURCE LINES 139-140 The error still contains the estimate of the error. .. GENERATED FROM PYTHON SOURCE LINES 140-142 .. code-block:: Python error[0] .. rst-class:: sphx-glr-script-out .. code-block:: none 1.6447801423976603e-09 .. GENERATED FROM PYTHON SOURCE LINES 143-150 During the algorithm, a collection of subintegrals .. math:: \int_{a_i}^{b_i} f(x) dx are approximated. The outputs :math:`a_i` and :math:`b_i` contain the subintervals used in the algorithm. .. GENERATED FROM PYTHON SOURCE LINES 150-153 .. code-block:: Python print("ai:", ai) print("bi:", bi) .. rst-class:: sphx-glr-script-out .. code-block:: none ai: [0,0.5,0.25,0.75,0.375,0.125,0.625,0.875,0.5625,0.9375,0.1875,0.3125]#12 bi: [0.125,0.5625,0.3125,0.875,0.5,0.1875,0.75,0.9375,0.625,1,0.25,0.375]#12 .. GENERATED FROM PYTHON SOURCE LINES 154-156 The corresponding value of the integrals are in :math:`f_i`. Since :math:`f` can be a multidimensional point, this is a :class:`~openturns.Sample`, which dimension corresponds to the output dimension of the function :math:`f`. .. GENERATED FROM PYTHON SOURCE LINES 156-158 .. code-block:: Python print("fi:", fi) .. rst-class:: sphx-glr-script-out .. code-block:: none fi: 0 : [ 0.108212 ] 1 : [ 0.10137 ] 2 : [ 0.0523839 ] 3 : [ 0.132726 ] 4 : [ 0.108212 ] 5 : [ 0.108834 ] 6 : [ 0.132726 ] 7 : [ 0.0738243 ] 8 : [ 0.0738243 ] 9 : [ 0.10137 ] 10 : [ 0.0523839 ] 11 : [ 0.108834 ] .. GENERATED FROM PYTHON SOURCE LINES 159-160 The sum of these sub-integrals is the value of the integral: .. GENERATED FROM PYTHON SOURCE LINES 160-162 .. code-block:: Python sum(fi)[0] .. rst-class:: sphx-glr-script-out .. code-block:: none 1.1547005383792528 .. GENERATED FROM PYTHON SOURCE LINES 163-164 The estimated error of each integral is in :math:`e_i`: .. GENERATED FROM PYTHON SOURCE LINES 164-172 .. code-block:: Python number_of_intervals = ai.getDimension() print("number of intervals:", number_of_intervals) for i in range(number_of_intervals): print( f"Integral #{i} : [{ai[i]:.4f}, {bi[i]:.4f}] is {fi[i, 0]:.4f} with error {ei[i]:.3e}" ) otv.View.ShowAll() .. rst-class:: sphx-glr-script-out .. code-block:: none number of intervals: 12 Integral #0 : [0.0000, 0.1250] is 0.1082 with error 4.326e-13 Integral #1 : [0.5000, 0.5625] is 0.1014 with error 2.991e-13 Integral #2 : [0.2500, 0.3125] is 0.0524 with error 1.249e-16 Integral #3 : [0.7500, 0.8750] is 0.1327 with error 1.163e-09 Integral #4 : [0.3750, 0.5000] is 0.1082 with error 4.326e-13 Integral #5 : [0.1250, 0.1875] is 0.1088 with error 3.485e-12 Integral #6 : [0.6250, 0.7500] is 0.1327 with error 1.163e-09 Integral #7 : [0.8750, 0.9375] is 0.0738 with error 2.665e-15 Integral #8 : [0.5625, 0.6250] is 0.0738 with error 2.637e-15 Integral #9 : [0.9375, 1.0000] is 0.1014 with error 2.991e-13 Integral #10 : [0.1875, 0.2500] is 0.0524 with error 1.318e-16 Integral #11 : [0.3125, 0.3750] is 0.1088 with error 3.485e-12 .. _sphx_glr_download_auto_numerical_methods_general_methods_plot_estimate_integral_Gauss_Kronrod.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_estimate_integral_Gauss_Kronrod.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_estimate_integral_Gauss_Kronrod.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_estimate_integral_Gauss_Kronrod.zip `