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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
|
# Numerical Integration
## Why Numerical?
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. In engineering, numerical integration provides a systematic approach to approximate the integral of a function over a finite interval.
---
## 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.
[Add diagrams]
### Midpoint Rule
The midpoint rule divides the interval into $n$ sub-intervals of equal width $h = (b-a)/n$ and
evaluates the function at the midpoint of each sub-interval:
$$
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)$.
```python
import numpy as np
a = 0
b = np.pi
n = 11
h = (b - a) / (n - 1)
x = np.linspace(a, b, n)
f = np.sin(x)
I_trap = (h/2)*(f[0] + 2 * sum(f[1:n-1]) + f[n-1])
err_trap = 2 - I_trap
print(I_trap)
print(err_trap)
```
### 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 sub-intervals $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].
$$
```python
import numpy as np
a = 0
b = np.pi
n = 11
h = (b - a) / (n - 1)
x = np.linspace(a, b, n)
f = np.sin(x)
I_simp = (h/3) * (f[0] + 2*sum(f[:n-2:2]) \
+ 4*sum(f[1:n-1:2]) + f[n-1])
err_simp = 2 - I_simp
print(I_simp)
print(err_simp)
```
**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.
```python
import numpy as np
a = 0
b = np.pi
n = 11
h = (b - a) / (n - 1)
x = np.linspace(a, b, n)
f = np.sin(x)
I_simp = (h/3) * (f[0] + 2*sum(f[:n-2:2]) \
+ 4*sum(f[1:n-1:2]) + f[n-1])
err_simp = 2 - I_simp
print(I_simp)
print(err_simp)
```
### 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.
```python
import numpy as np
a = 0.0
b = 1.0
f = lambda x: np.exp(-x**2)
iter_ = 0
iter_max = 50
err = 1.0
tol = 1e-14
I = np.zeros((iter_max + 2, iter_max + 2), dtype=float)
def trap_func(func, a, b, n):
"""Composite trapezoidal rule with n subintervals on [a, b]."""
h = (b - a) / n
x = np.linspace(a, b, n + 1)
y = func(x)
return h * (0.5 * y[0] + y[1:-1].sum() + 0.5 * y[-1])
n = 1
I[0, 0] = trap_func(f, a, b, n)
# ---- Romberg loop ----
while (iter_ < iter_max) and (err > tol):
iter_ += 1
n = 2 ** iter_
I[iter_, 0] = trap_func(f, a, b, n)
for k in range(1, iter_ + 1):
j = iter_ - k
I[j, k] = (4**k * I[j + 1, k - 1] - I[j, k - 1]) / (4**k - 1)
if iter_ >= 1:
err = abs((I[0, iter_] - I_]()
```
### Gaussian Quadrature
Gaussian quadrature chooses evaluation points and weights optimally 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.
```python
import numpy as np
a, b = 0.0, 1.0 # integration limits [a, b]
n_list = [2, 3, 4, 5, 8, 16] # quadrature orders to try
for n in n_list:
xi, wi = np.polynomial.legendre.leggauss(n)
x = 0.5*(b - a)*xi + 0.5*(b + a)
f_vals = np.exp(-x**2) # f(x) = e^(-x^2)
I_approx = 0.5*(b - a) * np.sum(wi * f_vals)
print(f"{n:<3}| {I_approx:.12f}")
```
```python
import numpy as np
from scipy.integrate import quad # for reference "exact" integral
# ---- Establish a, b, and the function ----
a = 0.0
b = 1.0
f = lambda x: np.exp(-x**2) # example: f(x) = e^(-x^2)
# ---- Common factor ----
dx = (b - a) / 2.0
# ---- Two-point Gauss-Legendre ----
x0 = ((b + a) + (b - a) * (-1/np.sqrt(3))) / 2
x1 = ((b + a) + (b - a) * ( 1/np.sqrt(3))) / 2
I_2Point = (f(x0) + f(x1)) * dx
print(f"2-point Gauss-Legendre: {I_2Point:.10f}")
# ---- Three-point Gauss-Legendre ----
c0, c1, c2 = 5/9, 8/9, 5/9
x0 = ((b + a) + (b - a) * (-np.sqrt(3/5))) / 2
x1 = ((b + a) + (b - a) * (0.0)) / 2
x2 = ((b + a) + (b - a) * ( np.sqrt(3/5))) / 2
I_3Point = (c0*f(x0) + c1*f(x1) + c2*f(x2)) * dx
print(f"3-point Gauss-Legendre: {I_3Point:.10f}")
# ---- Actual value using scipy.integrate.quad ----
I_actual, _ = quad(f, a, b)
print(f"Actual (scipy.quad): {I_actual:.10f}")
# ---- Errors ----
err2 = abs(I_actual - I_2Point)
err3 = abs(I_actual - I_3Point)
print(f"Error (2-point): {err2:.2e}")
print(f"Error (3-point): {err3:.2e}")
```
[[3_4_More_integral]]
|
