<jupyter_start><jupyter_text>Integration<jupyter_code>import numpy as np

from integration_algorithms import monte_carlo_quasi_two_dimensions
from integration_algorithms import quadrature_newton_trapezoid_one
from integration_algorithms import quadrature_newton_simpson_one
from integration_algorithms import quadrature_gauss_legendre_one
from integration_algorithms import quadrature_gauss_legendre_two
from integration_algorithms import monte_carlo_naive_one

from integration_plots import plot_naive_monte_carlo_randomness
from integration_plots import plot_naive_monte_carlo_error
from integration_plots import plot_gauss_legendre_weights
from integration_plots import plot_benchmarking_exercise
from integration_plots import plot_naive_monte_carlo
from integration_plots import plot_quasi_monte_carlo
from integration_plots import plot_starting_illustration
from integration_plots import plot_trapezoid_rule_illustration
from integration_plots import plot_simpsons_rule_illustration

from integration_problems import problem_kinked
from integration_problems import problem_smooth<jupyter_output><empty_output><jupyter_text>Outline1. Setup2. Newton-Cotes rules3. Gaussian formulas4. Monte Carlo integration5. Resources  SetupConsider finding the area under a continuous real-valued function $f$ over a bounded interval $[a, b]$:$$\begin{aligned}A = \int_a^b f(x) dx\end{aligned}$$Let's look a visual representation of our problem.<jupyter_code>plot_starting_illustration()<jupyter_output><empty_output><jupyter_text>Most numerical methods for computing this integral split up the original integral into a sum of several integrals, each covering a smaller part of the original integration interval $[a, b]$. This re-writing of the integral is based on a selection of integration points $x_i , i = 0, 1, . . . , n$ that are distributed on the interval $[a, b]$. Integration points may, or may not, be evenly distributed.Given the integration points, the original integral is re-written as a sum of integrals, each integral being computed over the sub-interval between two consecutive integration points. The integral from the beginning is thus expressed as:$$\begin{aligned}\int_a^b f(x) dx = \int_{x_0}^{x_1} f(x) dx + \int_{x_1}^{x_2} f(x) dx +  ...  + \int_{x_{n-1}}^{x_n} f(x) dx , \end{aligned}$$where $x_0 = a$ and $x_n = b$.The different integration methods will differ in the way they approximate each integral on the right hand side. The fundamental idea is that each term is an integral over a small interval $[x_i , x_{i+1}]$, and over this small interval, it makes sense to approximate $f$ by a simple shape  Newton-Cotes rulesNewton–Cotes quadrature rules are a group of formulas for numerical integration (also called quadrature) based on evaluating the integrand at equally spaced points. The integration points are then computed as$$\begin{aligned}x i = a + ih, i = 0, 1, . . . , n,\end{aligned}$$where $h = (b - a) / n$The closed Newton-Cotes formula of degree $n$ is stated as$$\begin{aligned}\int_a^b f(x) dx \approx \sum^n_{i=0} \omega_i f(x_i).\end{aligned}$$We will consider the first two degrees in more detail:* Trapezoid rule* Simpson’s rule  Trapezoid ruleThe Trapezoid rule approximates the area under the function $f$ with the area under a piecewise linear approximation to $f$.<jupyter_code>plot_trapezoid_rule_illustration()<jupyter_output><empty_output><jupyter_text>More formally:$$\begin{aligned}\int_a^b f(x) dx & = \int_{x_0}^{x_1} f(x) dx + \int_{x_1}^{x_2} f(x) dx  +  ... + \int_{x_{n-1}}^{x_n} f(x) dx \\ &\approx h \frac{f(x_0) - f(x_1)}{2} +  h \frac{f(x_1) - f(x_2)}{2} + ... + h \frac{f(x_{n-1}) - f(x_n)}{2}, \end{aligned}$$which we can further simplify to:$$\begin{aligned}\int_a^b f(x) dx \approx 0.5h f(x_0) + h\sum_{i = 1}^{n -1}f(x_i) + 0.5 h f(x_n).\end{aligned}$$The weights are then the following:$$\begin{aligned}\omega_i = \begin{cases}h / 2 & \text{if $i \in \{0, n\}$}   \\h & \text{otherwise}\end{cases}\end{aligned}$$<jupyter_code>??quadrature_newton_trapezoid_one
integrand = quadrature_newton_trapezoid_one(np.exp, 0, 1, 1000)
np.testing.assert_almost_equal(integrand, np.exp(1) - 1)<jupyter_output><empty_output><jupyter_text>Simpson’s ruleSimpson’s rule approximates the area under a function $f$ with the area under a piecewise quadratic approximation to $f$.<jupyter_code>plot_simpsons_rule_illustration()<jupyter_output><empty_output><jupyter_text>The weights are the following:$$\begin{aligned}w_i  = \begin{cases}h / 3 & \text{for $i \in [0, n]$} \\4 h / 3& \text{for $0 < i < n, i$ even} \\2 h / 3& \text{for $0 < i < n, i$ odd} . \\\end{cases}\end{aligned}$$<jupyter_code>??quadrature_newton_simpson_one
integrand = quadrature_newton_simpson_one(np.exp, 0, 1, 1001)
np.testing.assert_almost_equal(integrand, np.exp(1) - 1)<jupyter_output><empty_output><jupyter_text>Gaussian formulasGaussian quadrature formulas employ a very different logic to compute the area under a curve $f$ over a bounded interval $[a, b]$. Specifically, the $n$ quadrature nodes $x_i$ and $n$ quadrature weights $\omega_i$ are chosen to exactly integrate polynomials of degree $2n - 1$ or less.Gaussian quadrature approximations of the form:$$\begin{aligned}\int_a^b f(x) w(x) dx \approx \sum_{i = 1}^{n}\omega _i f(x_i)\end{aligned}$$for any nonnegative weighing function $w(x)$. Gauss-Legendre uses the weighing function $w(x) = 1$. Other versions of $w(x)$ lead to include Gauss-Chebyshev, Gauss-Hermite and Gauss-Laguerre.We start from the special domain between $[-1, 1]$, again the formula for our approximate solution looks very similar.$$\begin{aligned}\int_{-1}^1 f(x) dx\approx \sum^n_{i=1} w_i f(x_i)\end{aligned}$$So, more generally using a change of variables:$$\begin{aligned}\int_a^b f(x) = \approx \frac{b -a }{2} \sum^n_{i=1} w_i f\left(\frac{(x_i + 1)(b - a)}{2} + a\right)\end{aligned}$$Unlike in the case of the Newton-Cotes approach, Gaussian  nodes are not uniformly spaced and do not include the integration limits.<jupyter_code>[plot_gauss_legendre_weights(deg) for deg in [5, 11, 21]]<jupyter_output><empty_output><jupyter_text>_Question_* What are the properties of the weights?<jupyter_code>??quadrature_gauss_legendre_one
integrand = quadrature_gauss_legendre_one(np.exp, 0, 1, 1000)
np.testing.assert_almost_equal(integrand, np.exp(1) - 1)<jupyter_output><empty_output><jupyter_text>Benchmarking _Exercise_* Compare the performance of our integration routines on the following two integrals.We study the following two cases:$$\begin{aligned}f(x) & = \int_{-1}^1 e^{-x} dx \\f(x) & = \int_{-1}^1 \sqrt{|x|} dx .\end{aligned}$$<jupyter_code>plot_benchmarking_exercise()<jupyter_output><empty_output><jupyter_text>How can we extend the ideas so far to multidimensional integration? $$\begin{aligned}\int_{x}\int_{y} f(x, y) dx dy .\end{aligned}$$The product rule approximates the integral with the following sum:$$\begin{aligned}\sum_{i_x = 1}^m \sum_{i_y = 1}^m \omega_{i_x} \omega_{i_y} f(x_{1} y_{1}) .\end{aligned}$$A difficulty of this approach is the curse of dimensionality as the number of nodes inreases exponentially with the number of dimensions.<jupyter_code>??quadrature_gauss_legendre_two<jupyter_output><empty_output><jupyter_text>Monte Carlo integrationNaive Monte Carlo integration uses random sampling of a function to numerically compute an estimate of its integral. We can approximate this integral by averaging $N$ samples of the function $f$ drawn from a uniform distribution between $a$ and $b$.$$\begin{aligned}\int_a^b f(x) dx \approx (b - a) \frac{1}{N} \sum_{i = 1}^N f(x_i)\end{aligned}$$Why?$$\begin{aligned}E\left[ (b - a)  \frac{1}{N} \sum_{i = 1}^N f(x_i)\right] & = (b - a)  \frac{1}{N}  \sum_{i = 1}^N E\left[f(x_i)\right] \\& =  (b - a)  \frac{1}{N}  \sum_{i = 1}^N\left[\int_a^b f(x) \frac{1}{b - a} dx\right] \\& =  \frac{1}{N}  \sum_{i = 1}^N \int_a^b f(x) dx  \\& = \int_a^b f(x) dx.  \end{aligned}$$<jupyter_code>??monte_carlo_naive_one<jupyter_output><empty_output><jupyter_text>Monte Carlo integration only converges very slowly. Let's consider our smooth benchmarking example from earlier and compare the naive Monte Carlo approach to the quadrature methods.<jupyter_code>plot_naive_monte_carlo_error(max_nodes=200)<jupyter_output><empty_output><jupyter_text>Naive Monte Carlo integration differs from our earlier methods as the approximation itself is a random variable. How about the level of randomness at a given number $(n = 50)$ of integration nodes?<jupyter_code>plot_naive_monte_carlo_randomness()<jupyter_output><empty_output><jupyter_text>However, Monte Carlo integration is particular useful when tackling multidimensional integrals. Let's consider the computation of a two dimensional integral going forward.<jupyter_code>plot_naive_monte_carlo(100)<jupyter_output><empty_output><jupyter_text>Quasi-Monte Carlo employ deterministic sequences of nodes $x_i$ with the property that$$\begin{aligned}\lim_{n \rightarrow \infty} \frac{b - a}{n} \sum_{i=1}^n f(x_i)= \int_{a}^b f(x) dx.\end{aligned}$$Deterministic sequences of nodes chosen to fill space in a regular manner typically provide more accurate integration approximations than pseudo-random sequences. These sequences are readily available in the [chaospy](chaospy.readthedocs.io) package. The documentation of the package also includes a tutorial on [Monte Carlo integration](https://chaospy.readthedocs.io/en/master/tutorials/monte_carlo_integration.html).<jupyter_code>plot_quasi_monte_carlo(3)<jupyter_output><empty_output><jupyter_text>Otherwise, all looks pretty similar to our earlier implementation. Note that the default rule actually takes us back to the naive Monte Carlo method.<jupyter_code>??monte_carlo_quasi_two_dimensions<jupyter_output><empty_output><jupyter_text>Let's consider the following test function from Genz (1984).   _Exercises_1. What are the properties of the function?2. Compare the performance naive and quasi Monte Carlo integration routines for $u = (0.5, 0.5)$ and $a = (5, 5)$ over the two-dimensional unit cube.3. How does Gauss-Legendre quadrature perform? Applying the integration routines to functions with other input arguments can be achieved with some small modifications. Here we pass in our benchmarking function without any default arguments to a modified naive Monte Carlo integration routine.<jupyter_code>def genz_discontinuous_no_defaults(x, u, a):
    if x[0] > u[0] or x[1] > u[1]:
        return 0
    else:
        return np.exp((a * x).sum())


def monte_carlo_naive_with_args(f, args=(), a=0, b=1, n=10, seed=128):
    np.random.seed(seed)
    xvals = np.random.uniform(low=a, high=b, size=2 * n).reshape(n, 2)
    volume = (b - a) ** 2

    fvals = np.tile(np.nan, n)
    weights = np.tile(1 / n, n)

    for i, xval in enumerate(xvals):
        # Here is the main modification that passes in the arguments by position now.
        fvals[i] = f(xval, *args)

    return volume * np.sum(weights * fvals)


u, a = (0.5, 0.5), (5, 5)
rslt = monte_carlo_naive_with_args(genz_discontinuous_no_defaults, args=(u, a))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11562415380449499, 'fraction_numerical': 0.01380991064175467, 'mean_word_length': 5.7975460122699385, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 71, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '291984', 'license': 'MIT', 'n_tokens_mistral': 3824, 'n_tokens_neox': 3581, 'n_words': 1120, 'path': 'lectures/integration/notebook.ipynb', 'repo_name': 'CasualDan/ose-course-scientific-computing', 'size': '11079'}