

# Basic numerical integration methods: Riemann sums and trapezoidal rule

Consider a function $y(x)$ defined on the interval $[a, b]$. The goal of numerical integration is to approximate the definite integral $\int_a^b y(x) dx$.

One basic numerical integration method is the Riemann sum, which approximates the area under the curve by dividing the interval into a large number of small subintervals and summing the areas of the rectangles under the curve. 

The trapezoidal rule is an improvement over the Riemann sum, which approximates the area under the curve by dividing the interval into a large number of small subintervals and summing the areas of the trapezoids under the curve.

The trapezoidal rule formula is:

$$\int_a^b y(x) dx \approx \frac{h}{2} \sum_{i=1}^n (y(x_i) + y(x_{i-1}))$$

where $h = \frac{b-a}{n}$ is the width of each subinterval, $x_i = a + ih$, and $n$ is the number of subintervals.

Consider the function $y(x) = x^2$ on the interval $[0, 1]$. We want to approximate the integral $\int_0^1 x^2 dx$ using the trapezoidal rule with $n = 5$.

Using the trapezoidal rule formula, we can approximate the integral as:

$$\int_0^1 x^2 dx \approx \frac{1}{2} \sum_{i=1}^5 (x_i^2 + x_{i-1}^2)$$

where $h = \frac{1}{5}$ and $x_i = \frac{i}{5}$.

## Exercise
Calculate the integral $\int_0^1 x^2 dx$ using the trapezoidal rule with $n = 10$.

The trapezoidal rule has a higher order of accuracy than the Riemann sum, but it still has some limitations. For example, it assumes that the function is continuous and that the intervals are of equal width.

# Adaptive integration techniques

Adaptive integration techniques are a class of numerical integration methods that adjust the width of the subintervals based on the error estimate. This allows them to achieve higher accuracy with fewer function evaluations.

One example of an adaptive integration technique is the adaptive Simpson's rule, which combines the trapezoidal rule and Simpson's rule. It starts with the trapezoidal rule and refines the intervals where the error is large.

The adaptive Simpson's rule formula is:

$$\int_a^b y(x) dx \approx \frac{h}{3} \sum_{i=1}^n (y(x_i) + 4y(x_{i-1}) + y(x_{i-2}))$$

where $h = \frac{b-a}{n}$ is the width of each subinterval, $x_i = a + ih$, and $n$ is the number of subintervals.

Consider the function $y(x) = x^2$ on the interval $[0, 1]$. We want to approximate the integral $\int_0^1 x^2 dx$ using the adaptive Simpson's rule with $n = 5$.

Using the adaptive Simpson's rule formula, we can approximate the integral as:

$$\int_0^1 x^2 dx \approx \frac{1}{3} \sum_{i=1}^5 (x_i^2 + 4x_{i-1}^2 + x_{i-2}^2)$$

where $h = \frac{1}{5}$ and $x_i = \frac{i}{5}$.

## Exercise
Calculate the integral $\int_0^1 x^2 dx$ using the adaptive Simpson's rule with $n = 10$.

Adaptive integration techniques are more accurate than fixed-step methods like the trapezoidal rule, but they require more function evaluations. This trade-off can be controlled by adjusting the error tolerance and the maximum number of function evaluations.

# Error estimation in numerical integration

To estimate the error in numerical integration, we can use the Richardson extrapolation method. This method refines the estimate by applying Richardson's extrapolation formula to the results of the integration method at different orders of accuracy.

The Richardson extrapolation formula is:

$$E_{n+1} = \frac{1}{2} (E_n - E_{n-1})$$

where $E_n$ is the estimate at order $n$ and $E_{n-1}$ is the estimate at order $n-1$.

Consider the function $y(x) = x^2$ on the interval $[0, 1]$. We want to approximate the integral $\int_0^1 x^2 dx$ using the trapezoidal rule with $n = 5$ and estimate the error using Richardson extrapolation.

Using the trapezoidal rule formula, we can approximate the integral as:

$$\int_0^1 x^2 dx \approx \frac{1}{2} \sum_{i=1}^5 (y(x_i) + y(x_{i-1}))$$

where $h = \frac{1}{5}$ and $x_i = \frac{i}{5}$.

Next, we can estimate the error using Richardson extrapolation.

## Exercise
Calculate the integral $\int_0^1 x^2 dx$ using the trapezoidal rule with $n = 5$ and estimate the error using Richardson extrapolation.

Error estimation is an important aspect of numerical integration, as it allows us to assess the accuracy of our results and make informed decisions about the number of function evaluations and the width of the subintervals.

# Quadrature methods: Simpson's rule, Gauss-Kronrod rule, and other advanced techniques

Quadrature methods are a class of numerical integration methods that use specialized formulas to approximate the integral. These methods can achieve high accuracy with fewer function evaluations compared to adaptive integration techniques.

One example of a quadrature method is Simpson's rule, which combines the trapezoidal rule and the midpoint rule. It approximates the integral by summing the areas of the parabolic segments under the curve.

The Simpson's rule formula is:

$$\int_a^b y(x) dx \approx \frac{h}{3} \sum_{i=1}^n (y(x_i) + 4y(x_{i-1}) + y(x_{i-2}))$$

where $h = \frac{b-a}{n}$ is the width of each subinterval, $x_i = a + ih$, and $n$ is the number of subintervals.

Consider the function $y(x) = x^2$ on the interval $[0, 1]$. We want to approximate the integral $\int_0^1 x^2 dx$ using Simpson's rule with $n = 5$.

Using the Simpson's rule formula, we can approximate the integral as:

$$\int_0^1 x^2 dx \approx \frac{1}{3} \sum_{i=1}^5 (x_i^2 + 4x_{i-1}^2 + x_{i-2}^2)$$

where $h = \frac{1}{5}$ and $x_i = \frac{i}{5}$.

## Exercise
Calculate the integral $\int_0^1 x^2 dx$ using Simpson's rule with $n = 10$.

Another example of a quadrature method is the Gauss-Kronrod rule, which combines the Gaussian quadrature and the Kronrod extension. It approximates the integral by summing the areas of the Gaussian quadrature segments under the curve and the Kronrod extension segments under the curve.

# Using SciPy's integration functions

SciPy is a popular Python library that provides a wide range of mathematical functions, including integration functions. These functions can be used to perform numerical integration on functions defined in Python.

One example of a SciPy integration function is `scipy.integrate.quad`, which approximates a definite integral using a user-specified integration method.

Consider the function $y(x) = x^2$ on the interval $[0, 1]$. We want to approximate the integral $\int_0^1 x^2 dx$ using the trapezoidal rule with `scipy.integrate.quad`.

```python
import scipy.integrate as spi

def y(x):
    return x**2

result, error = spi.quad(y, 0, 1, method='trapz')

print(f'Integral: {result}, Error: {error}')
```

## Exercise
Calculate the integral $\int_0^1 x^2 dx$ using the trapezoidal rule with `scipy.integrate.quad`.

SciPy's integration functions provide a convenient and powerful way to perform numerical integration in Python. They can be used to integrate functions defined in Python or C, and they support a wide range of integration methods and options.

# Applications of numerical integration in the SciPy package

SciPy's integration functions can be used for a wide range of applications in the SciPy package. Some examples include:

- Solving ordinary differential equations with numerical integration.
- Integrating functions with complex domains.
- Handling discontinuities and singularities in functions.
- Comparison of numerical integration methods.

# Solving ordinary differential equations with numerical integration

Numerical integration can be used to solve ordinary differential equations (ODEs) by approximating the solution as a function of the independent variable.

One example of a numerical integration method for solving ODEs is the Euler method, which approximates the solution as a weighted average of the function values at consecutive points.

Consider the first-order ordinary differential equation $y'(x) = f(x, y)$, where $y(x_0) = y_0$. We want to approximate the solution using the Euler method with step size $h$.

```python
def f(x, y):
    return x + y

x0, y0 = 0, 1
x = np.arange(x0, x0 + 10, h)

y = np.zeros(len(x))
y[0] = y0

for i in range(1, len(x)):
    y[i] = y[i-1] + h * f(x[i-1], y[i-1])

print(y)
```

## Exercise
Solve the first-order ordinary differential equation $y'(x) = x + y$ with initial condition $y(0) = 1$ using the Euler method with step size $h = 0.1$.

Numerical integration can be used to approximate the solutions of various types of ordinary differential equations, including nonlinear equations and systems of equations.

# Integrating functions with complex domains

Numerical integration can also be used to approximate the integral of functions defined on complex domains. This can be useful in various applications, such as solving partial differential equations, simulating physical systems, and analyzing complex data.

One example of a numerical integration method for complex domains is the complex trapezoidal rule, which approximates the integral by summing the areas of the complex trapezoids under the curve.

Consider the function $y(x) = x^2$ on the complex interval $[0, 1]$ in the complex plane. We want to approximate the integral $\int_0^1 x^2 dx$ using the complex trapezoidal rule with $n = 5$.

Using the complex trapezoidal rule formula, we can approximate the integral as:

$$\int_0^1 x^2 dx \approx \frac{1}{2} \sum_{i=1}^5 (y(x_i) + y(x_{i-1}))$$

where $h = \frac{1}{5}$ and $x_i = \frac{i}{5}$.

## Exercise
Calculate the integral $\int_0^1 x^2 dx$ using the complex trapezoidal rule with $n = 10$.

Numerical integration with complex domains can be used to approximate the integral of functions defined on various complex domains, including non-rectangular domains and domains with singularities.

# Handling discontinuities and singularities in functions

Numerical integration can be used to approximate the integral of functions with discontinuities and singularities. This can be challenging, as the accuracy of the numerical integration methods may be affected by the presence of these features.

One approach to handling discontinuities and singularities is to use adaptive integration techniques that can adaptively refine the intervals where the function is discontinuous or singular.

Consider the function $y(x) = \frac{1}{x}$ on the interval $[0, 1]$. We want to approximate the integral $\int_0^1 \frac{1}{x} dx$ using the adaptive Simpson's rule with $n = 5$.

Using the adaptive Simpson's rule formula, we can approximate the integral as:

$$\int_0^1 \frac{1}{x} dx \approx \frac{1}{3} \sum_{i=1}^5 (y(x_i) + 4y(x_{i-1}) + y(x_{i-2}))$$

where $h = \frac{1}{5}$ and $x_i = \frac{i}{5}$.

## Exercise
Calculate the integral $\int_0^1 \frac{1}{x} dx$ using the adaptive Simpson's rule with $n = 10$.

Handling discontinuities and singularities in functions can be challenging, but advanced numerical integration methods like adaptive integration techniques can help overcome these challenges and provide accurate results.

# Comparison of numerical integration methods

There are many numerical integration methods available in SciPy and other libraries. The choice of the method depends on the specific requirements of the problem, such as the accuracy, the number of function evaluations, and the complexity of the function.

Some common numerical integration methods include:

- Trapezoidal rule.
- Adaptive Simpson's rule.
- Simpson's rule.
- Gauss-Kronrod rule.
- Euler method.
- Complex trapezoidal rule.

Each method has its advantages and disadvantages, and the choice of the method should be based on the specific needs of the problem.

## Exercise
Compare the accuracy and efficiency of the trapezoidal rule, adaptive Simpson's rule, and Simpson's rule for approximating the integral of the function $y(x) = x^2$ on the interval $[0, 1]$.

Comparing numerical integration methods is an important aspect of numerical integration, as it allows us to choose the most suitable method for a given problem and ensure the accuracy and efficiency of the results.