---
author:
- "by [Eugeniy E. Mikhailov](http://physics.wm.edu/~evmik/) and [Greg
  Bentsen](https://physics.wm.edu/~gbentsen/)"
title: Numerical integration
---

## Logistics and Agenda

## Midterm 1:

Do it at home part

-   Midterm write up and code is due Sunday October 4th 11:59pm

-   You may use AI freely to write code, debug, and test.

-   **You are responsible for the code you submit.**

-   We expect you to understand what the code is doing.

-   You may **not** use AI to write the report.

    -   Report should be done by you and only one per group is needed

In class part

-   on Monday October 5th, with computers but no AI and no Internet
    connection
    -   we will give you modified assignment which would be similar
    -   it will require modification or reuse of prepared code

# Integration problem statement

Suppose we are given function $$f(x)$$ our goal is to find
$$\int_a^b f(x) dx$$ Not all function can be easily integrated
analytically in the elementary enough form.

::: example
$$\int_0^y \frac{\exp(-x^2)}{x^3+\cos(\omega x)} dx$$
:::

So we must use numerical methods.

## The Riemann sum

Recall the Riemann integral definition our goal is to find

$$\int_a^b f(x) dx = \lim_{N \to \infty} \sum_{i=1}^{N-1} f(x_i) h$$

where $N$ is the number of points, $h=(b-a)/(N-1)$ is the distance
between points.

![Function under integral](./pics/func2integrate.png){width="70%"}

## The Rectangle (box) method

Riemann rule is almost direct recipe for the rectangle left end method

$$\int_a^b f(x) dx \approx \sum_{i=1}^{N-1} f(x_i) h$$ where
$h=\frac{b-a}{N-1}$ and $x_i=a+(i-1)h$

We just need to remember about round off errors so $h$ should not be too
small or equivalently $N$ should not be to big.

![Rectangle left end method](./pics/rectangle.png){width="70%"}

## Trapezoidal method

$$\int_a^b f(x) dx \approx h \times (\frac{1}{2} f_1 + f_2 +f_3 + \cdots
    +f_{N-2} +f_{N-1} + \frac{1}{2} f_{N}) = h \sum_{i=1}^N f(x_i) w_i$$

where $h=\frac{b-a}{N-1}$ and $x_i=a+(i-1)h$

![Trapezoid method](./pics/trapezoidal.png){width="80%"}

## Let's compare implementation of trapezoidal method

::: columns
::: {.column width="50% "}
-   machine generated with prompt "use python and write integration by
    trapezoid method, it should take a, b, and Np - number of points"

see file [trapezoid_ai.py](./code/trapezoid_ai.py)

``` phython

# AI generated code

def trapezoidal_ai(func, a, b, Np):
    """
    Numerical integration using the trapezoidal rule.

    Parameters:
    func: function to integrate
    a: lower limit of integration
    b: upper limit of integration
    Np: number of points (intervals = Np - 1)

    Returns:
    Approximation of the integral
    """

# Calculate step size

h = (b - a) / (Np - 1)

# Evaluate function at endpoints

integral = 0.5 * (func(a) + func(b))

# Sum the interior points

for i in range(1, Np - 1):
        x = a + i * h
        integral += func(x)

# Multiply by step size

integral *= h

    return integral
```
:::

::: {.column width="50%"}
-   human programmed

see file [trapezoid_human.py](./code/trapezoid_human.py)

``` phython

# human generated code

import numpy as np

def trapezoidal_human(func, a, b, Np):
    """Numerical integration using the trapezoidal rule."""
    h = (b - a) / (Np - 1)

    x = np.linspace(a, b, Np)
    y = func(x)
    y[0] = y[0] / 2
    y[-1] = y[-1] / 2

    return np.sum(y) * h
```

This code also runs `\alert{faster by factor of 30}`{=tex}
:::
:::

## Simpson method: approximation by parabolas.

Main idea take two adjacent points $f_{k-1}$, $f_{k}$, and $f_{k+1}$ and
approximate the $f(x)$ by parabola on this $2h$ interval. Then the
integral on this interval is approximated as
$\frac{h}{3} (f_{k-1} + 4f_{k} + f_{k+1})$

$$\int_a^b f(x) dx \approx h \frac{1}{3}\times ( f_1 + 4 f_2 + 2 f_3 + 4
        f_4 + \cdots +2 f_{N-2} +4 f_{N-1} + f_{N})
        =  h \sum_{i=1}^N f(x_i) w_i $$

where $h=\frac{b-a}{N-1}$ and $x_i=a+(i-1)h$

Note that N must be in special form N=2k+1, i.e. odd.

![Simpson method](./pics/simpson.png){width="70%"}

## Integration error estimate

Rectangle method
$$E={\mathbf{O}}\left( \frac{ (b -a) h}{2} f'  \right) = 
{\mathbf{O}}\left( \frac{ (b -a)^2}{2 N} f'  \right)$$

Trapezoidal method
$$E={\mathbf{O}}\left( \frac{ (b -a) h^2}{12} f''  \right) = 
            {\mathbf{O}}\left( \frac{ (b -a)^3}{12 N^2} f''  \right)$$

Simpson method
$$E={\mathbf{O}}\left( \frac{ (b -a) h^4}{180} f^{(4)}  \right) = 
{\mathbf{O}}\left( \frac{ (b -a)^5}{180 N^4} f^{(4)}  \right)$$

::: example
Here $\mathbf{O}$ represents [big O
notation](https://en.wikipedia.org/wiki/Big_O_notation) for so called
"order of approximation".

$$ f(x) = \mathbf{O(x)}$$ means that $|f(x)|  \le C|x|$, where $C$ is a
constant
:::

## Library method for numerical integration

use `quad` from `scipy.integrate` module

``` python
import scipy.integrate as integrate

def f(x): return x*x
int_estimate, error_estimate = integrate.quad(f, 0, 1)

In [2]: int_estimate
Out[2]: 0.33333333333333337

In [3]: error_estimate
Out[3]: 3.700743415417189e-15
```