.. 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_iterated_quadrature.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_iterated_quadrature.py: Estimate a multivariate integral with IteratedQuadrature ======================================================== .. GENERATED FROM PYTHON SOURCE LINES 7-19 Introduction ------------ In this example, we consider a function :math:`f: \Rset^d \mapsto \Rset` and we compute the following integral: .. math:: I_f = \int_{a}^{b}\, \int_{l_1(x_0)}^{u_1(x_0)}\, \int_{l_2(x_0, x_1)}^{u_2(x_0,x_1)}\, \int_{l_{d-1}(x_0, \dots, x_{d-2})}^{u_{d-1}(x_0, \dots, x_{d-2})} \, f(x_0, \dots, x_{d-1})\mathrm{d}{x_{d-1}}\dots\mathrm{d}{x_0} using an iterated quadrature algorithm :class:`~openturns.IteratedQuadrature`, based on the Gauss-Kronrod :class:`~openturns.GaussKronrod` 1d quadrature algorithm. .. GENERATED FROM PYTHON SOURCE LINES 20-25 .. code-block:: Python import openturns as ot import openturns.viewer as otv import math as m .. GENERATED FROM PYTHON SOURCE LINES 26-27 To get better colours for the function iso-lines. .. GENERATED FROM PYTHON SOURCE LINES 27-29 .. code-block:: Python ot.ResourceMap.SetAsString("Contour-DefaultColorMapNorm", "rank") .. GENERATED FROM PYTHON SOURCE LINES 30-45 Case 1: The integrand function presents a lot of discontinuities. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ We consider the function: .. math:: f : \begin{array}{lcl} \Rset^2 & \rightarrow & \Rset \\ (x,y) & \rightarrow & 1_{{x^2+y^2 \leq 4}}(x,y) \end{array} the domain :math:`\set{D} =[-2,2]\times [-2,2]` and the integral: .. math:: I = \int_{\set{D}} f(x,y)\,dxdy = 4\pi .. GENERATED FROM PYTHON SOURCE LINES 47-48 We first define the integrand :math:`f` and the domain :math:`\set{D}`. .. GENERATED FROM PYTHON SOURCE LINES 48-60 .. code-block:: Python def kernel1(x): if x[0] ** 2 + x[1] ** 2 <= 4: return [1.0] else: return [0.0] integrand = ot.PythonFunction(2, 1, kernel1) domain = ot.Interval([-2.0, -2.0], [2.0, 2.0]) .. GENERATED FROM PYTHON SOURCE LINES 61-63 We draw the iso-lines of the integrand function which is constant equal to 1 inside the circle with radius equal to 2 and equal to 0.0 outside. .. GENERATED FROM PYTHON SOURCE LINES 63-69 .. code-block:: Python g = integrand.draw([-3.0, -3.0], [3.0, 3.0]) g.setTitle(r"Function $f(x,y) = 1_{\{x^2+y^2 \leq 4\}}(x,y)$") g.setXTitle("x") g.setYTitle("y") view = otv.View(g) .. image-sg:: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_001.png :alt: Function $f(x,y) = 1_{\{x^2+y^2 \leq 4\}}(x,y)$ :srcset: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 70-80 We compute the integral using an :class:`~openturns.IteratedQuadrature` algorithm with the :class:`~openturns.GaussKronrodRule` *G11K23*, which means that we use on each subinterval: - 11 points for the Gauss approximations, - 23 points for the Kronrod approximation including the 11 previous ones. We limit the number of subintervals and we define a maximum error between both approximations: when the difference between both approximations is lower, a new subdivision is performed. .. GENERATED FROM PYTHON SOURCE LINES 80-86 .. code-block:: Python maximumSubIntervals = 32 maximumError = 1e-4 GKRule = ot.GaussKronrodRule(ot.GaussKronrodRule.G1K3) integ_algo = ot.GaussKronrod(maximumSubIntervals, maximumError, GKRule) integration_algo = ot.IteratedQuadrature(integ_algo) .. GENERATED FROM PYTHON SOURCE LINES 87-89 To get the nodes used to approximate the integral, we need to use a :class:`~openturns.MemoizeFunction` that stores all the calls to the function. .. GENERATED FROM PYTHON SOURCE LINES 89-96 .. code-block:: Python integrand_memoized = ot.MemoizeFunction(integrand) I_value = integration_algo.integrate(integrand_memoized, domain) print("I = ", I_value[0]) print("Exact value = ", 4 * m.pi) nodes = integrand_memoized.getInputHistory() print("Nodes : ", nodes) .. rst-class:: sphx-glr-script-out .. code-block:: none I = 12.898438724470733 Exact value = 12.566370614359172 Nodes : 0 : [ -1 -1 ] 1 : [ -1 -1.7746 ] 2 : [ -1 -0.225403 ] ... 11337 : [ 1.68574 1.5 ] 11338 : [ 1.68574 1.1127 ] 11339 : [ 1.68574 1.8873 ] .. GENERATED FROM PYTHON SOURCE LINES 97-100 We can draw superimpose the nodes on the graph with the iso-lines of the function :math:`f`. We can see that the algorithm focuses on the nodes where the function has its discontinuities. .. GENERATED FROM PYTHON SOURCE LINES 100-108 .. code-block:: Python cloud_nodes = ot.Cloud(nodes) cloud_nodes.setPointStyle("dot") cloud_nodes.setColor("blue") g.add(cloud_nodes) g.setLegendPosition("") view = otv.View(g) .. image-sg:: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_002.png :alt: Function $f(x,y) = 1_{\{x^2+y^2 \leq 4\}}(x,y)$ :srcset: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 109-129 Case 2: The integrand function has large gradients. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ We consider the function: .. math:: f : \begin{array}{lcl} \Rset^2 & \rightarrow & \Rset \\ (x,y) & \rightarrow & e^{-(x^2+y^2)} + e^{-8((x-4)^2+(y-4)^2)} \end{array} the domain :math:`\set{D} =[-2,2]\times [-2,2]` and the integral: .. math:: I = \int_{\set{D}} f(x,y)\, dxdy. Using Maple, we obtain the reference value: .. math:: I_{ref} = 3.51961338289448753812591427600 .. GENERATED FROM PYTHON SOURCE LINES 131-132 We first define the integrand :math:`f` and the domain :math:`\set{D}`. .. GENERATED FROM PYTHON SOURCE LINES 132-144 .. code-block:: Python def kernel2(x): return [ m.exp(-(x[0] ** 2 + x[1] ** 2)) + m.exp(-8 * ((x[0] - 4) ** 2 + (x[1] - 4) ** 2)) ] integrand = ot.PythonFunction(2, 1, kernel2) domain = ot.Interval([-2.0, -2.0], [6.0, 6.0]) .. GENERATED FROM PYTHON SOURCE LINES 145-146 We draw the iso-lines of the integrand function. .. GENERATED FROM PYTHON SOURCE LINES 146-151 .. code-block:: Python g = integrand.draw([-3.0, -3.0], [7.0, 7.0]) g.setTitle(r"Function $f(x,y) = e^{-(x^2+y^2)} + e^{-8((x-4)^2+(y-4)^2)} $") g.setXTitle("x") g.setYTitle("y") .. GENERATED FROM PYTHON SOURCE LINES 152-154 To get the nodes used to approximate the integral, we need to use a :class:`~openturns.MemoizeFunction` that stores all the calls to the function. .. GENERATED FROM PYTHON SOURCE LINES 154-160 .. code-block:: Python integrand_memoized = ot.MemoizeFunction(integrand) I_value = integration_algo.integrate(integrand_memoized, domain) print("I = ", I_value[0]) nodes = integrand_memoized.getInputHistory() print("Nodes : ", nodes) .. rst-class:: sphx-glr-script-out .. code-block:: none I = 3.5196013466116494 Nodes : 0 : [ 0 0 ] 1 : [ 0 -1.54919 ] 2 : [ 0 1.54919 ] ... 30165 : [ 3.48591 4.04688 ] 30166 : [ 3.48591 4.03477 ] 30167 : [ 3.48591 4.05898 ] .. GENERATED FROM PYTHON SOURCE LINES 161-164 We can draw in blue the nodes on the graph that contains the iso-lines of the function :math:`f`. We can see that the algorithm focuses on the nodes where the function has fast variations. .. GENERATED FROM PYTHON SOURCE LINES 164-173 .. code-block:: Python cloud_nodes = ot.Cloud(nodes) cloud_nodes.setPointStyle("dot") cloud_nodes.setColor("blue") g.add(cloud_nodes) g.setLegendPosition("") view = otv.View(g) .. image-sg:: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_003.png :alt: Function $f(x,y) = e^{-(x^2+y^2)} + e^{-8((x-4)^2+(y-4)^2)} $ :srcset: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 174-195 Case 3: The integration domain is complex. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ We consider the function: .. math:: f : \begin{array}{lcl} \Rset^2 & \rightarrow & \Rset \\ (x,y) & \rightarrow & \cos(2x) \sin(1.5y) \end{array} the domain :math:`\set{D} = \{(x,y)\, |\, -2\pi \leq x \leq 2\pi, -2 - \cos(x) \leq y \leq 2 + \cos(x) \}` and the integral: .. math:: I = \int_{\set{D}} f(x,y)\, dxdy. Using Maple, we obtain the reference value: .. math:: I_{ref} = -0.548768615494004063851543284777 .. GENERATED FROM PYTHON SOURCE LINES 197-198 We first define the integrand :math:`f` and the domain :math:`\set{D}`. .. GENERATED FROM PYTHON SOURCE LINES 198-206 .. code-block:: Python def kernel3(x): return [m.cos(2.0 * x[0]) * m.cos(1.5 * x[1])] integrand = ot.PythonFunction(2, 1, kernel3) .. GENERATED FROM PYTHON SOURCE LINES 207-217 The integration domain is defined by the two following functions: .. math:: x \rightarrow u(x) = 2 + \cos(x) .. math:: x \rightarrow \ell(x) = -2 - \cos(x) .. GENERATED FROM PYTHON SOURCE LINES 217-220 .. code-block:: Python upper_func = ot.SymbolicFunction(["x"], [" 2 + cos(x)"]) lower_func = ot.SymbolicFunction(["x"], ["-2 - cos(x)"]) .. GENERATED FROM PYTHON SOURCE LINES 221-224 We draw the iso-lines of the integrand function in the domain of integration. To do that, we define a new function which is the restriction of the integrand function to the integration domain. .. GENERATED FROM PYTHON SOURCE LINES 224-242 .. code-block:: Python def kernel3_insideDomain(x): low_x = lower_func([x[0]])[0] up_x = upper_func([x[0]])[0] if x[1] > low_x and x[1] < up_x: return kernel3(x) else: return [0.0] integrand_domain = ot.PythonFunction(2, 1, kernel3_insideDomain) a = 2 * m.pi b = 4.0 # Here, we increase the default number of levels. ot.ResourceMap.SetAsUnsignedInteger("Contour-DefaultLevelsNumber", 30) g = integrand_domain.draw([-a, -b], [a, b]) .. GENERATED FROM PYTHON SOURCE LINES 243-244 We add the bounds of the domain :math:`\set{D}`. .. GENERATED FROM PYTHON SOURCE LINES 244-259 .. code-block:: Python low_bound = lower_func.draw([-a, -b], [a, b]).getDrawable(0) up_bound = upper_func.draw([-a, -b], [a, b]).getDrawable(0) low_bound.setColor("red") low_bound.setLineWidth(2) up_bound.setColor("red") up_bound.setLineWidth(2) g.add(low_bound) g.add(up_bound) g.setTitle(r"Function $f(x,y) = cos(2x)sin(1.5y)$") g.setXTitle("x") g.setYTitle("y") view = otv.View(g) # sphinx_gallery_thumbnail_number = 4 .. image-sg:: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_004.png :alt: Function $f(x,y) = cos(2x)sin(1.5y)$ :srcset: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 260-262 To get the nodes used to approximate the integral, we need to use a :class:`~openturns.MemoizeFunction` that stores all the calls to the function. .. GENERATED FROM PYTHON SOURCE LINES 262-276 .. code-block:: Python integrand_memoized = ot.MemoizeFunction(integrand) maximumSubIntervals = 4 maximumError = 1e-4 GKRule = ot.GaussKronrodRule(ot.GaussKronrodRule.G7K15) integration_algo = ot.IteratedQuadrature( ot.GaussKronrod(maximumSubIntervals, maximumError, GKRule) ) I_value = integration_algo.integrate( integrand_memoized, -a, a, [lower_func], [upper_func] ) print("I = ", I_value[0]) nodes = integrand_memoized.getInputHistory() print("Nodes : ", nodes) .. rst-class:: sphx-glr-script-out .. code-block:: none I = -0.5487686157122131 Nodes : 0 : [ -3.14159 -0.5 ] 1 : [ -3.14159 -0.995728 ] 2 : [ -3.14159 -0.00427231 ] ... 2697 : [ 5.03878 1.63122 ] 2698 : [ 5.03878 0.919216 ] 2699 : [ 5.03878 1.40141 ] .. GENERATED FROM PYTHON SOURCE LINES 277-280 We can superimpose the nodes on the graph with the iso-lines of the function :math:`f`. We can see that the algorithm focuses on the nodes where the function has fast variations. .. GENERATED FROM PYTHON SOURCE LINES 280-288 .. code-block:: Python cloud_nodes = ot.Cloud(nodes) cloud_nodes.setPointStyle("dot") cloud_nodes.setColor("black") g.add(cloud_nodes) g.setLegendPosition("") view = otv.View(g) .. image-sg:: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_005.png :alt: Function $f(x,y) = cos(2x)sin(1.5y)$ :srcset: /auto_numerical_methods/general_methods/images/sphx_glr_plot_estimate_integral_iterated_quadrature_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 289-290 Reset ResourceMap .. GENERATED FROM PYTHON SOURCE LINES 290-292 .. code-block:: Python ot.ResourceMap.Reload() .. GENERATED FROM PYTHON SOURCE LINES 293-294 Show all the graphs. .. GENERATED FROM PYTHON SOURCE LINES 294-295 .. code-block:: Python view.ShowAll() .. _sphx_glr_download_auto_numerical_methods_general_methods_plot_estimate_integral_iterated_quadrature.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_iterated_quadrature.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_estimate_integral_iterated_quadrature.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_estimate_integral_iterated_quadrature.zip `