1"""
2์์น ์ ๋ถ (Numerical Integration)
3Numerical Integration Methods
4
5์ ์ ๋ถ โซ[a,b] f(x)dx ๋ฅผ ์์น์ ์ผ๋ก ๊ณ์ฐํ๋ ๋ฐฉ๋ฒ๋ค์
๋๋ค.
6"""
7
8import numpy as np
9import matplotlib.pyplot as plt
10from typing import Callable, Tuple
11
12
13# =============================================================================
14# 1. ์ง์ฌ๊ฐํ ๋ฒ (Rectangular Rule)
15# =============================================================================
16def rectangular_left(f: Callable, a: float, b: float, n: int) -> float:
17 """
18 ์ผ์ชฝ ์ง์ฌ๊ฐํ ๋ฒ
19 ๊ฐ ๊ตฌ๊ฐ์ ์ผ์ชฝ ๋์ ์์ ํจ์๊ฐ์ผ๋ก ์ง์ฌ๊ฐํ ๋์ด ๊ฒฐ์
20 """
21 h = (b - a) / n
22 result = 0
23 for i in range(n):
24 x = a + i * h
25 result += f(x)
26 return h * result
27
28
29def rectangular_right(f: Callable, a: float, b: float, n: int) -> float:
30 """์ค๋ฅธ์ชฝ ์ง์ฌ๊ฐํ ๋ฒ"""
31 h = (b - a) / n
32 result = 0
33 for i in range(1, n + 1):
34 x = a + i * h
35 result += f(x)
36 return h * result
37
38
39def rectangular_midpoint(f: Callable, a: float, b: float, n: int) -> float:
40 """
41 ์ค์ ์ง์ฌ๊ฐํ ๋ฒ (Midpoint Rule)
42 ์ค์ฐจ: O(hยฒ)
43 """
44 h = (b - a) / n
45 result = 0
46 for i in range(n):
47 x = a + (i + 0.5) * h
48 result += f(x)
49 return h * result
50
51
52# =============================================================================
53# 2. ์ฌ๋ค๋ฆฌ๊ผด ๋ฒ (Trapezoidal Rule)
54# =============================================================================
55def trapezoidal(f: Callable, a: float, b: float, n: int) -> float:
56 """
57 ์ฌ๋ค๋ฆฌ๊ผด ๋ฒ
58 ๊ฐ ๊ตฌ๊ฐ์ ์ฌ๋ค๋ฆฌ๊ผด๋ก ๊ทผ์ฌ
59 ์ค์ฐจ: O(hยฒ)
60
61 ๊ณต์: h/2 * [f(xโ) + 2*ฮฃf(xแตข) + f(xโ)]
62 """
63 h = (b - a) / n
64 result = f(a) + f(b)
65
66 for i in range(1, n):
67 x = a + i * h
68 result += 2 * f(x)
69
70 return h * result / 2
71
72
73def trapezoidal_vectorized(f: Callable, a: float, b: float, n: int) -> float:
74 """์ฌ๋ค๋ฆฌ๊ผด ๋ฒ (NumPy ๋ฒกํฐํ)"""
75 x = np.linspace(a, b, n + 1)
76 y = f(x)
77 return np.trapz(y, x)
78
79
80# =============================================================================
81# 3. ์ฌํ์จ ๋ฒ (Simpson's Rule)
82# =============================================================================
83def simpsons_rule(f: Callable, a: float, b: float, n: int) -> float:
84 """
85 ์ฌํ์จ 1/3 ๋ฒ์น
86 ๊ฐ ๊ตฌ๊ฐ์ 2์ฐจ ๋คํญ์์ผ๋ก ๊ทผ์ฌ
87 ์ค์ฐจ: O(hโด) - ๋งค์ฐ ์ ํ
88
89 ์กฐ๊ฑด: n์ ์ง์์ฌ์ผ ํจ
90 ๊ณต์: h/3 * [f(xโ) + 4*ฮฃf(ํ์) + 2*ฮฃf(์ง์) + f(xโ)]
91 """
92 if n % 2 != 0:
93 n += 1 # ์ง์๋ก ์กฐ์
94
95 h = (b - a) / n
96 result = f(a) + f(b)
97
98 for i in range(1, n):
99 x = a + i * h
100 if i % 2 == 0:
101 result += 2 * f(x)
102 else:
103 result += 4 * f(x)
104
105 return h * result / 3
106
107
108def simpsons_38(f: Callable, a: float, b: float, n: int) -> float:
109 """
110 ์ฌํ์จ 3/8 ๋ฒ์น
111 3์ฐจ ๋คํญ์ ๊ทผ์ฌ
112 ์กฐ๊ฑด: n์ 3์ ๋ฐฐ์์ฌ์ผ ํจ
113 """
114 if n % 3 != 0:
115 n = (n // 3 + 1) * 3
116
117 h = (b - a) / n
118 result = f(a) + f(b)
119
120 for i in range(1, n):
121 x = a + i * h
122 if i % 3 == 0:
123 result += 2 * f(x)
124 else:
125 result += 3 * f(x)
126
127 return 3 * h * result / 8
128
129
130# =============================================================================
131# 4. ๋กฌ๋ฒ ๋ฅด๊ทธ ์ ๋ถ (Romberg Integration)
132# =============================================================================
133def romberg(f: Callable, a: float, b: float, max_order: int = 5) -> float:
134 """
135 ๋กฌ๋ฒ ๋ฅด๊ทธ ์ ๋ถ
136 ์ฌ๋ค๋ฆฌ๊ผด ๋ฒ์ Richardson ์ธ์ฝ์ ๋ฐ๋ณต ์ ์ฉ
137 ๋งค์ฐ ๋์ ์ ํ๋ ๋ฌ์ฑ ๊ฐ๋ฅ
138 """
139 R = np.zeros((max_order, max_order))
140
141 # ์ฒซ ๋ฒ์งธ ์ด: ์ฌ๋ค๋ฆฌ๊ผด ๋ฒ
142 h = b - a
143 R[0, 0] = h * (f(a) + f(b)) / 2
144
145 for i in range(1, max_order):
146 h = h / 2
147 # ์๋ก์ด ์ ๋ค์ ํฉ
148 sum_new = sum(f(a + (2*k - 1) * h) for k in range(1, 2**(i-1) + 1))
149 R[i, 0] = R[i-1, 0] / 2 + h * sum_new
150
151 # Richardson ์ธ์ฝ
152 for j in range(1, i + 1):
153 R[i, j] = R[i, j-1] + (R[i, j-1] - R[i-1, j-1]) / (4**j - 1)
154
155 return R[max_order-1, max_order-1]
156
157
158# =============================================================================
159# 5. ๊ฐ์ฐ์ค ๊ตฌ์ ๋ฒ (Gaussian Quadrature)
160# =============================================================================
161def gauss_legendre(f: Callable, a: float, b: float, n: int = 5) -> float:
162 """
163 ๊ฐ์ฐ์ค-๋ฅด์ฅ๋๋ฅด ๊ตฌ์ ๋ฒ
164 n๊ฐ์ ์ ์ผ๋ก (2n-1)์ฐจ ๋คํญ์๊น์ง ์ ํํ ์ ๋ถ
165 ๋งค์ฐ ํจ์จ์
166
167 [-1, 1]์์ [a, b]๋ก ๋ณํ ์ ์ฉ
168 """
169 # n๊ฐ ๋
ธ๋์ ๊ฐ์ค์น (์ฌ์ ๊ณ์ฐ๋ ๊ฐ)
170 nodes, weights = np.polynomial.legendre.leggauss(n)
171
172 # ๊ตฌ๊ฐ ๋ณํ: [-1, 1] โ [a, b]
173 # x = (b-a)/2 * t + (a+b)/2
174 # dx = (b-a)/2 * dt
175
176 transformed_nodes = (b - a) / 2 * nodes + (a + b) / 2
177 result = sum(w * f(x) for x, w in zip(transformed_nodes, weights))
178
179 return (b - a) / 2 * result
180
181
182# =============================================================================
183# 6. ์ ์์ ์ ๋ถ (Adaptive Quadrature)
184# =============================================================================
185def adaptive_simpsons(
186 f: Callable,
187 a: float,
188 b: float,
189 tol: float = 1e-8,
190 max_depth: int = 50
191) -> Tuple[float, int]:
192 """
193 ์ ์์ ์ฌํ์จ ์ ๋ถ
194 ์ค์ฐจ๊ฐ ํฐ ๊ตฌ๊ฐ์ ์ฌ๊ท์ ์ผ๋ก ์ธ๋ถํ
195 """
196 call_count = [0]
197
198 def _adaptive(a, b, fa, fb, fc, S, tol, depth):
199 call_count[0] += 1
200 c = (a + b) / 2
201 d = (a + c) / 2
202 e = (c + b) / 2
203
204 fd = f(d)
205 fe = f(e)
206
207 S_left = (c - a) / 6 * (fa + 4*fd + fc)
208 S_right = (b - c) / 6 * (fc + 4*fe + fb)
209 S_new = S_left + S_right
210
211 # ์ค์ฐจ ์ถ์
212 error = (S_new - S) / 15
213
214 if depth >= max_depth or abs(error) <= tol:
215 return S_new + error # Richardson ์ธ์ฝ
216 else:
217 left = _adaptive(a, c, fa, fc, fd, S_left, tol/2, depth+1)
218 right = _adaptive(c, b, fc, fb, fe, S_right, tol/2, depth+1)
219 return left + right
220
221 fa, fb = f(a), f(b)
222 fc = f((a + b) / 2)
223 S = (b - a) / 6 * (fa + 4*fc + fb)
224
225 result = _adaptive(a, b, fa, fb, fc, S, tol, 0)
226 return result, call_count[0]
227
228
229# =============================================================================
230# ์ค์ฐจ ๋ถ์ ๋ฐ ์๊ฐํ
231# =============================================================================
232def analyze_convergence(f: Callable, a: float, b: float, exact: float):
233 """์๋ ด ์๋ ๋ถ์"""
234 ns = [4, 8, 16, 32, 64, 128, 256, 512]
235 methods = {
236 'Midpoint': rectangular_midpoint,
237 'Trapezoidal': trapezoidal,
238 'Simpson': simpsons_rule,
239 }
240
241 print("\n์๋ ด ๋ถ์:")
242 print("-" * 70)
243 print(f"{'n':>6} | {'Midpoint':>14} | {'Trapezoidal':>14} | {'Simpson':>14}")
244 print("-" * 70)
245
246 errors = {name: [] for name in methods}
247
248 for n in ns:
249 row = f"{n:>6} |"
250 for name, method in methods.items():
251 result = method(f, a, b, n)
252 error = abs(result - exact)
253 errors[name].append(error)
254 row += f" {error:>14.2e} |"
255 print(row)
256
257 return ns, errors
258
259
260def plot_methods_comparison(f, a, b, n=10):
261 """์ ๋ถ ๋ฐฉ๋ฒ ์๊ฐํ"""
262 fig, axes = plt.subplots(2, 2, figsize=(12, 10))
263
264 x_dense = np.linspace(a, b, 500)
265 y_dense = f(x_dense)
266
267 methods = [
268 ('์ง์ฌ๊ฐํ (์ค์ )', rectangular_midpoint),
269 ('์ฌ๋ค๋ฆฌ๊ผด', trapezoidal),
270 ('์ฌํ์จ', simpsons_rule),
271 ]
272
273 # ํจ์์ ์ ๋ถ ์์ญ
274 for ax, (name, method) in zip(axes.flat[:3], methods):
275 ax.plot(x_dense, y_dense, 'b-', linewidth=2, label='f(x)')
276 ax.fill_between(x_dense, y_dense, alpha=0.3)
277 ax.axhline(y=0, color='k', linewidth=0.5)
278
279 h = (b - a) / n
280 x_pts = np.linspace(a, b, n + 1)
281
282 if 'mid' in name:
283 for i in range(n):
284 xm = a + (i + 0.5) * h
285 ax.bar(xm, f(xm), width=h, alpha=0.5, edgecolor='r', fill=False)
286 elif 'trap' in name.lower() or '์ฌ๋ค๋ฆฌ๊ผด' in name:
287 for i in range(n):
288 x0, x1 = x_pts[i], x_pts[i+1]
289 ax.fill([x0, x1, x1, x0], [0, 0, f(x1), f(x0)], alpha=0.5, edgecolor='r', fill=False)
290
291 result = method(f, a, b, n)
292 ax.set_title(f'{name}: {result:.6f}')
293 ax.set_xlabel('x')
294 ax.set_ylabel('y')
295 ax.grid(True, alpha=0.3)
296
297 # ์๋ ด ๊ทธ๋ํ
298 ax = axes[1, 1]
299 ns, errors = [4, 8, 16, 32, 64, 128], {name: [] for name, _ in methods}
300 exact = 2.0 # โซ[0,ฯ] sin(x)dx = 2
301
302 for n in ns:
303 for name, method in methods:
304 errors[name].append(abs(method(f, a, b, n) - exact))
305
306 for name, errs in errors.items():
307 ax.loglog(ns, errs, 'o-', label=name)
308
309 ax.set_xlabel('n (๊ตฌ๊ฐ ์)')
310 ax.set_ylabel('์ค์ฐจ')
311 ax.set_title('์๋ ด ์๋')
312 ax.legend()
313 ax.grid(True, alpha=0.3)
314
315 plt.tight_layout()
316 plt.savefig('/opt/projects/01_Personal/03_Study/Numerical_Simulation/examples/numerical_integration.png', dpi=150)
317 plt.close()
318 print("๊ทธ๋ํ ์ ์ฅ: numerical_integration.png")
319
320
321# =============================================================================
322# ํ
์คํธ
323# =============================================================================
324def main():
325 print("=" * 60)
326 print("์์น ์ ๋ถ (Numerical Integration) ์์ ")
327 print("=" * 60)
328
329 # ์์ 1: โซ[0,ฯ] sin(x)dx = 2
330 print("\n[์์ 1] โซ[0,ฯ] sin(x)dx = 2")
331 print("-" * 50)
332
333 f = np.sin
334 a, b = 0, np.pi
335 exact = 2.0
336
337 n = 100
338 results = {
339 '์ค์ ์ง์ฌ๊ฐํ': rectangular_midpoint(f, a, b, n),
340 '์ฌ๋ค๋ฆฌ๊ผด': trapezoidal(f, a, b, n),
341 '์ฌํ์จ 1/3': simpsons_rule(f, a, b, n),
342 '๋กฌ๋ฒ ๋ฅด๊ทธ': romberg(f, a, b, 6),
343 '๊ฐ์ฐ์ค-๋ฅด์ฅ๋๋ฅด (5์ )': gauss_legendre(f, a, b, 5),
344 }
345
346 print(f"์ ํ๊ฐ: {exact}")
347 print(f"{'๋ฐฉ๋ฒ':<20} {'๊ฒฐ๊ณผ':<15} {'์ค์ฐจ':<15}")
348 print("-" * 50)
349 for name, result in results.items():
350 error = abs(result - exact)
351 print(f"{name:<20} {result:<15.10f} {error:<15.2e}")
352
353 # ์์ 2: โซ[0,1] e^(-xยฒ)dx โ 0.7468...
354 print("\n[์์ 2] โซ[0,1] e^(-xยฒ)dx (๊ฐ์ฐ์ค ์ ๋ถ)")
355 print("-" * 50)
356
357 f2 = lambda x: np.exp(-x**2)
358 exact2 = 0.746824132812427 # ์ฐธ๊ฐ (erf ์ฌ์ฉ)
359
360 for n in [10, 50, 100]:
361 trap = trapezoidal(f2, 0, 1, n)
362 simp = simpsons_rule(f2, 0, 1, n)
363 print(f"n={n:3d}: ์ฌ๋ค๋ฆฌ๊ผด={trap:.10f}, ์ฌํ์จ={simp:.10f}")
364
365 gauss = gauss_legendre(f2, 0, 1, 10)
366 print(f"๊ฐ์ฐ์ค-๋ฅด์ฅ๋๋ฅด (10์ ): {gauss:.10f}")
367 print(f"์ ํ๊ฐ: {exact2:.10f}")
368
369 # ์์ 3: ์ ์์ ์ ๋ถ
370 print("\n[์์ 3] ์ ์์ ์ ๋ถ (๊ธ๋ณํ๋ ํจ์)")
371 print("-" * 50)
372
373 f3 = lambda x: 1 / (1 + 100 * (x - 0.5)**2) # ์ข์ ํผํฌ
374 exact3 = 0.3141277802329 # โ ฯ/10
375
376 trap = trapezoidal(f3, 0, 1, 100)
377 result, calls = adaptive_simpsons(f3, 0, 1, tol=1e-8)
378
379 print(f"์ฌ๋ค๋ฆฌ๊ผด (n=100): {trap:.10f}, ์ค์ฐจ: {abs(trap - exact3):.2e}")
380 print(f"์ ์์ ์ฌํ์จ: {result:.10f}, ์ค์ฐจ: {abs(result - exact3):.2e}, ํธ์ถ: {calls}")
381
382 # ์๊ฐํ
383 try:
384 plot_methods_comparison(np.sin, 0, np.pi, 10)
385 except Exception as e:
386 print(f"๊ทธ๋ํ ์์ฑ ์คํจ: {e}")
387
388 # ์๋ ด ๋ถ์
389 print("\n" + "=" * 60)
390 print("์๋ ด ์๋ ๋ถ์ (โซ[0,ฯ] sin(x)dx)")
391 analyze_convergence(np.sin, 0, np.pi, 2.0)
392
393 print("\n" + "=" * 60)
394 print("์์น ์ ๋ถ ๋ฐฉ๋ฒ ๋น๊ต")
395 print("=" * 60)
396 print("""
397 | ๋ฐฉ๋ฒ | ์ค์ฐจ ์ฐจ์ | ํน์ง |
398 |--------------|----------|-------------------------------|
399 | ์ง์ฌ๊ฐํ (์ข/์ฐ)| O(h) | ๊ฐ์ฅ ๋จ์, ๋ถ์ ํ |
400 | ์ค์ ์ง์ฌ๊ฐํ | O(hยฒ) | ์ง์ฌ๊ฐํ ์ค ๊ฐ์ฅ ์ ํ |
401 | ์ฌ๋ค๋ฆฌ๊ผด | O(hยฒ) | ๋จ์ํ๊ณ ํจ์จ์ |
402 | ์ฌํ์จ 1/3 | O(hโด) | ๋งค์ฐ ์ ํ, ๋ถ๋๋ฌ์ด ํจ์์ ์ ํฉ |
403 | ๋กฌ๋ฒ ๋ฅด๊ทธ | ~O(h^2k) | Richardson ์ธ์ฝ, ๊ณ ์ ํ๋ |
404 | ๊ฐ์ฐ์ค ๊ตฌ์ | ~O(h^2n) | ์ต์ ์ ์ผ๋ก ์ต๋ ์ ํ๋ |
405 | ์ ์์ ์ ๋ถ | ๊ฐ๋ณ | ๊ธ๋ณ ๊ตฌ๊ฐ ์๋ ์ธ๋ถํ |
406
407 ์ค๋ฌด ๊ถ์ฅ:
408 - scipy.integrate.quad: ์ ์์ ๊ฐ์ฐ์ค ๊ตฌ์ (๊ฐ์ฅ ์ผ๋ฐ์ )
409 - scipy.integrate.romberg: ๋กฌ๋ฒ ๋ฅด๊ทธ ์ ๋ถ
410 - scipy.integrate.simps: ์ฌํ์จ ๋ฒ์น (๋ฑ๊ฐ๊ฒฉ ๋ฐ์ดํฐ)
411 - scipy.integrate.trapz: ์ฌ๋ค๋ฆฌ๊ผด (๋ฑ๊ฐ๊ฒฉ ๋ฐ์ดํฐ)
412 """)
413
414
415if __name__ == "__main__":
416 main()