summaryrefslogtreecommitdiff
path: root/tutorials/module_3/numerical_integration.md
blob: 86ce015769eac521bf01562b3e93bd3bab4fe3a4 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
## Midpoint Method


## Trapezoidal Method


## Romberg Integration


## Gaussian Integration


## Simpson's Rule

### Simpsons 1/3

### Simpsons 3/8




# Numerical Integration
## Introduction and Importance
Integration is one of the fundamental tools in engineering analysis. Mechanical engineers frequently encounter integrals when computing work from force–displacement data, determining heat transfer from a time-dependent signal, or calculating lift and drag forces from pressure distributions over an airfoil. While some integrals can be evaluated analytically, most practical problems involve functions that are either too complex or are available only as experimental data. As engineering we choose numerical integration—also known as quadrature—provides a systematic approach to approximate the integral of a function over a finite interval.

In this tutorial, we will study several standard methods of numerical integration, compare their
accuracy, and implement them in Python. By the end, you will understand not only how to apply
each method, but also when one method may be more suitable than another.

---

## Section 2 — Core Concepts and Numerical Methods

### General Form of Numerical Integration
We wish to approximate a definite integral of the form:
$$
I = \int_a^b f(x) \, dx
$$
by a weighted sum of function values:
$$
I \approx \sum_{i=0}^m w_i f(x_i).
$$
Here, $x_i$ are the chosen evaluation points and $w_i$ are their associated weights.

### Midpoint Rule
The midpoint rule divides the interval into $n$ subintervals of equal width $h = (b-a)/n$ and
evaluates the function at the midpoint of each subinterval:
$$
I \approx \sum_{i=0}^{n-1} h \, f\!\left(x_i + \tfrac{h}{2}\right).
$$
This method achieves second-order accuracy (error decreases with \(h^2\)).

### Trapezoidal Rule
The trapezoidal rule approximates the area under the curve as a series of trapezoids:
$$
I \approx \frac{h}{2}\Big[f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n)\Big].
$$
It is simple to implement and works especially well for tabulated data. Like the midpoint rule,
its accuracy is of order \(O(h^2)\).

### Simpson’s Rule
Simpson’s rules use polynomial interpolation to achieve higher accuracy.

- **Simpson’s 1/3 Rule (order \(O(h^4)\))**
  Requires an even number of subintervals \(n\):
  $$
  I \approx \frac{h}{3}\Big[f(x_0) + 4\sum_{\text{odd } i} f(x_i) +
  2\sum_{\text{even } i<n} f(x_i) + f(x_n)\Big].
  $$

- **Simpson’s 3/8 Rule**
  Based on cubic interpolation across three subintervals:
  $$
  I \approx \frac{3h}{8}\Big[f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)\Big].
  $$
  Often used to handle the final three panels when \(n\) is not divisible by two.

### Romberg Integration
Romberg integration refines trapezoidal approximations using Richardson extrapolation to cancel
error terms. This results in very rapid convergence for smooth functions, making it a powerful
method once the trapezoidal rule is implemented.

### Gaussian Quadrature
Gaussian quadrature chooses evaluation points and weights optimally (based on Legendre
polynomials) to maximize accuracy. With \(n\) evaluation points, Gaussian quadrature is exact for
polynomials of degree up to \(2n-1\). This makes it extremely efficient for smooth integrands.

---

## Exercise — Implementing the Composite Trapezoidal Rule

**Objective:** Implement a Python function to approximate integrals using the trapezoidal rule.

```python
import numpy as np

def trapz(f, a, b, n):
    x = np.linspace(a, b, n+1)
    y = f(x)
    h = (b - a) / n
    return h * (0.5*y[0] + y[1:-1].sum() + 0.5*y[-1])

# Example tests
f1 = np.sin
I_true1 = 2.0   # ∫_0^π sin(x) dx
for n in [4, 8, 16, 32]:
    print(n, trapz(f1, 0, np.pi, n))
```

Students should compare results for increasing \(n\) and observe how the error decreases with
\(O(h^2)\).

---

## Section 4 — Extensions and Advanced Topics

- **Romberg Integration:** Demonstrates how extrapolation accelerates convergence of trapezoidal
approximations.
- **Gaussian Quadrature:** Introduces optimal integration points and highlights efficiency for
polynomial and smooth functions.
- **Simpson’s Rules:** Show how higher-order polynomial interpolation improves accuracy.

These methods will be implemented and compared in subsequent assignments to build a deeper
understanding of numerical integration accuracy and efficiency.

---

## Assignment 1 — Gaussian Quadrature

Implement a Python function for two-point and three-point Gauss–Legendre quadrature over an
arbitrary interval \([a,b]\). Verify exactness for polynomials up to the appropriate degree and
compare performance against the trapezoidal rule on oscillatory test functions.

---

## Assignment 2 — Simpson’s 1/3 Rule

Implement the composite Simpson’s 1/3 rule. Test its accuracy on smooth functions and compare its performance to the trapezoidal rule and Gaussian quadrature. Document error trends and discuss cases where Simpson’s method is preferable.

---