1"""
2λ£½κ²-μΏ ν λ°©λ² (Runge-Kutta Methods)
3Runge-Kutta Methods for ODEs
4
5κ³ μ°¨ μ νλλ₯Ό κ°μ§ ODE μμΉ ν΄λ²μ
λλ€.
6"""
7
8import numpy as np
9import matplotlib.pyplot as plt
10from typing import Callable, Tuple, List, Union
11
12
13# =============================================================================
14# 1. RK2 (2μ°¨ λ£½κ²-μΏ ν)
15# =============================================================================
16def rk2_midpoint(
17 f: Callable[[float, float], float],
18 y0: float,
19 t_span: Tuple[float, float],
20 h: float
21) -> Tuple[np.ndarray, np.ndarray]:
22 """
23 RK2 μ€μ λ°©λ² (Midpoint Method)
24
25 k1 = f(t_n, y_n)
26 k2 = f(t_n + h/2, y_n + h*k1/2)
27 y_{n+1} = y_n + h*k2
28
29 μ€μ°¨: O(hΒ²)
30 """
31 t0, tf = t_span
32 n_steps = int((tf - t0) / h)
33
34 t = np.linspace(t0, tf, n_steps + 1)
35 y = np.zeros(n_steps + 1)
36 y[0] = y0
37
38 for i in range(n_steps):
39 k1 = f(t[i], y[i])
40 k2 = f(t[i] + h/2, y[i] + h*k1/2)
41 y[i + 1] = y[i] + h * k2
42
43 return t, y
44
45
46def rk2_heun(
47 f: Callable[[float, float], float],
48 y0: float,
49 t_span: Tuple[float, float],
50 h: float
51) -> Tuple[np.ndarray, np.ndarray]:
52 """
53 RK2 Heun λ°©λ² (μμ μ€μΌλ¬μ λμΌ)
54
55 k1 = f(t_n, y_n)
56 k2 = f(t_n + h, y_n + h*k1)
57 y_{n+1} = y_n + h*(k1 + k2)/2
58
59 μ€μ°¨: O(hΒ²)
60 """
61 t0, tf = t_span
62 n_steps = int((tf - t0) / h)
63
64 t = np.linspace(t0, tf, n_steps + 1)
65 y = np.zeros(n_steps + 1)
66 y[0] = y0
67
68 for i in range(n_steps):
69 k1 = f(t[i], y[i])
70 k2 = f(t[i] + h, y[i] + h*k1)
71 y[i + 1] = y[i] + h * (k1 + k2) / 2
72
73 return t, y
74
75
76# =============================================================================
77# 2. RK4 (4μ°¨ λ£½κ²-μΏ ν) - κ°μ₯ λ리 μ¬μ©λ¨
78# =============================================================================
79def rk4(
80 f: Callable[[float, float], float],
81 y0: float,
82 t_span: Tuple[float, float],
83 h: float
84) -> Tuple[np.ndarray, np.ndarray]:
85 """
86 κ³ μ μ RK4 (Classical 4th-order Runge-Kutta)
87
88 k1 = f(t_n, y_n)
89 k2 = f(t_n + h/2, y_n + h*k1/2)
90 k3 = f(t_n + h/2, y_n + h*k2/2)
91 k4 = f(t_n + h, y_n + h*k3)
92 y_{n+1} = y_n + h*(k1 + 2*k2 + 2*k3 + k4)/6
93
94 μ€μ°¨: O(hβ΄)
95 κ°μ₯ λ리 μ¬μ©λλ λ°©λ²
96 """
97 t0, tf = t_span
98 n_steps = int((tf - t0) / h)
99
100 t = np.linspace(t0, tf, n_steps + 1)
101 y = np.zeros(n_steps + 1)
102 y[0] = y0
103
104 for i in range(n_steps):
105 k1 = f(t[i], y[i])
106 k2 = f(t[i] + h/2, y[i] + h*k1/2)
107 k3 = f(t[i] + h/2, y[i] + h*k2/2)
108 k4 = f(t[i] + h, y[i] + h*k3)
109 y[i + 1] = y[i] + h * (k1 + 2*k2 + 2*k3 + k4) / 6
110
111 return t, y
112
113
114# =============================================================================
115# 3. λ²‘ν° RK4 (μ°λ¦½ ODE)
116# =============================================================================
117def rk4_system(
118 f: Callable[[float, np.ndarray], np.ndarray],
119 y0: np.ndarray,
120 t_span: Tuple[float, float],
121 h: float
122) -> Tuple[np.ndarray, np.ndarray]:
123 """
124 μ°λ¦½ ODEλ₯Ό μν RK4
125
126 dy/dt = f(t, y), y = [y1, y2, ..., yn]
127 """
128 t0, tf = t_span
129 n_steps = int((tf - t0) / h)
130 n_vars = len(y0)
131
132 t = np.linspace(t0, tf, n_steps + 1)
133 y = np.zeros((n_steps + 1, n_vars))
134 y[0] = y0
135
136 for i in range(n_steps):
137 k1 = np.array(f(t[i], y[i]))
138 k2 = np.array(f(t[i] + h/2, y[i] + h*k1/2))
139 k3 = np.array(f(t[i] + h/2, y[i] + h*k2/2))
140 k4 = np.array(f(t[i] + h, y[i] + h*k3))
141 y[i + 1] = y[i] + h * (k1 + 2*k2 + 2*k3 + k4) / 6
142
143 return t, y
144
145
146# =============================================================================
147# 4. μ μμ RK45 (Runge-Kutta-Fehlberg)
148# =============================================================================
149def rkf45(
150 f: Callable[[float, float], float],
151 y0: float,
152 t_span: Tuple[float, float],
153 tol: float = 1e-6,
154 h_init: float = 0.1,
155 h_min: float = 1e-10,
156 h_max: float = 1.0
157) -> Tuple[np.ndarray, np.ndarray]:
158 """
159 μ μμ RK4(5) (Runge-Kutta-Fehlberg)
160
161 4μ°¨μ 5μ°¨ κ·Όμ¬λ₯Ό λμμ κ³μ°νμ¬ μ€μ°¨ μΆμ
162 μ€μ°¨μ λ°λΌ μ€ν
ν¬κΈ° μλ μ‘°μ
163 """
164 # RK45 κ³μ (Fehlberg)
165 c = [0, 1/4, 3/8, 12/13, 1, 1/2]
166 a = [
167 [],
168 [1/4],
169 [3/32, 9/32],
170 [1932/2197, -7200/2197, 7296/2197],
171 [439/216, -8, 3680/513, -845/4104],
172 [-8/27, 2, -3544/2565, 1859/4104, -11/40]
173 ]
174 b4 = [25/216, 0, 1408/2565, 2197/4104, -1/5, 0]
175 b5 = [16/135, 0, 6656/12825, 28561/56430, -9/50, 2/55]
176
177 t0, tf = t_span
178 t_list = [t0]
179 y_list = [y0]
180
181 t = t0
182 y = y0
183 h = h_init
184
185 while t < tf:
186 if t + h > tf:
187 h = tf - t
188
189 # 6κ°μ k κ³μ°
190 k = [0] * 6
191 k[0] = f(t, y)
192 for j in range(1, 6):
193 y_temp = y + h * sum(a[j][m] * k[m] for m in range(j))
194 k[j] = f(t + c[j] * h, y_temp)
195
196 # 4μ°¨μ 5μ°¨ κ·Όμ¬
197 y4 = y + h * sum(b4[j] * k[j] for j in range(6))
198 y5 = y + h * sum(b5[j] * k[j] for j in range(6))
199
200 # μ€μ°¨ μΆμ
201 error = abs(y5 - y4)
202
203 if error < 1e-15:
204 error = 1e-15
205
206 # μ€ν
ν¬κΈ° μ‘°μ
207 h_new = 0.9 * h * (tol / error) ** 0.2
208
209 if error <= tol:
210 # μ€ν
μ±κ³΅
211 t = t + h
212 y = y5
213 t_list.append(t)
214 y_list.append(y)
215 h = min(h_max, max(h_min, h_new))
216 else:
217 # μ€ν
μ€ν¨, μ¬μλ
218 h = max(h_min, h_new)
219
220 return np.array(t_list), np.array(y_list)
221
222
223# =============================================================================
224# 5. RK4 vs Euler λΉκ΅
225# =============================================================================
226def compare_methods():
227 """λ€μν λ°©λ² λΉκ΅"""
228 # dy/dt = y, y(0) = 1 β y = e^t
229 f = lambda t, y: y
230 y0 = 1
231 t_span = (0, 2)
232 exact = lambda t: np.exp(t)
233
234 print("\nRK λ°©λ² λΉκ΅ (dy/dt = y, t β [0, 2])")
235 print("-" * 70)
236 print(f"{'h':>10} | {'μ€μΌλ¬':>12} | {'RK2':>12} | {'RK4':>12}")
237 print("-" * 70)
238
239 from os.path import dirname, abspath
240 import sys
241 sys.path.insert(0, dirname(abspath(__file__)))
242
243 try:
244 # 03_ode_euler.pyλ μ«μλ‘ μμνμ¬ μ§μ import λΆκ°
245 from importlib import import_module
246 euler_module = import_module("03_ode_euler")
247 euler_forward = euler_module.euler_forward
248 except (ImportError, ModuleNotFoundError):
249 def euler_forward(f, y0, t_span, h):
250 t0, tf = t_span
251 n_steps = int((tf - t0) / h)
252 t = np.linspace(t0, tf, n_steps + 1)
253 y = np.zeros(n_steps + 1)
254 y[0] = y0
255 for i in range(n_steps):
256 y[i + 1] = y[i] + h * f(t[i], y[i])
257 return t, y
258
259 hs = [0.2, 0.1, 0.05, 0.025, 0.0125]
260
261 for h in hs:
262 _, y_euler = euler_forward(f, y0, t_span, h)
263 _, y_rk2 = rk2_heun(f, y0, t_span, h)
264 _, y_rk4 = rk4(f, y0, t_span, h)
265
266 e_euler = abs(y_euler[-1] - exact(2))
267 e_rk2 = abs(y_rk2[-1] - exact(2))
268 e_rk4 = abs(y_rk4[-1] - exact(2))
269
270 print(f"{h:>10.4f} | {e_euler:>12.2e} | {e_rk2:>12.2e} | {e_rk4:>12.2e}")
271
272
273# =============================================================================
274# μκ°ν
275# =============================================================================
276def plot_rk_comparison():
277 """RK λ°©λ² λΉκ΅ μκ°ν"""
278 fig, axes = plt.subplots(2, 2, figsize=(12, 10))
279
280 # μμ 1: μ§μ μ±μ₯
281 f = lambda t, y: y
282 exact = lambda t: np.exp(t)
283 t_exact = np.linspace(0, 2, 200)
284
285 ax = axes[0, 0]
286 ax.plot(t_exact, exact(t_exact), 'k-', linewidth=2, label='μ νν΄')
287 for h in [0.5, 0.25, 0.1]:
288 t, y = rk4(f, 1, (0, 2), h)
289 ax.plot(t, y, 'o-', markersize=4, label=f'RK4 h={h}')
290 ax.set_title("RK4: dy/dt = y")
291 ax.legend()
292 ax.grid(True)
293
294 # μμ 2: μ‘°ν μ§λμ (μλμ§ λ³΄μ‘΄)
295 def harmonic(t, state):
296 y, v = state
297 return np.array([v, -y])
298
299 ax = axes[0, 1]
300 t_exact = np.linspace(0, 20, 500)
301 ax.plot(t_exact, np.cos(t_exact), 'k-', linewidth=2, label='μ νν΄')
302
303 t, sol = rk4_system(harmonic, np.array([1, 0]), (0, 20), 0.1)
304 ax.plot(t, sol[:, 0], 'r-', label='RK4 h=0.1')
305
306 t, sol = rk4_system(harmonic, np.array([1, 0]), (0, 20), 0.01)
307 ax.plot(t, sol[:, 0], 'b--', label='RK4 h=0.01')
308
309 ax.set_title("μ‘°ν μ§λμ: y'' = -y")
310 ax.legend()
311 ax.grid(True)
312
313 # μμ 3: μ μμ vs κ³ μ μ€ν
314 f_stiff = lambda t, y: -50 * (y - np.cos(t))
315
316 ax = axes[1, 0]
317 t_ex = np.linspace(0, 1, 500)
318 y_ex = np.cos(t_ex) + (0 - 1) * np.exp(-50 * t_ex) # κ·Όμ¬ ν΄μν΄
319
320 t1, y1 = rk4(f_stiff, 0, (0, 1), 0.01)
321 t2, y2 = rkf45(f_stiff, 0, (0, 1), tol=1e-6)
322
323 ax.plot(t_ex, y_ex, 'k-', linewidth=2, label='μ°Έμ‘°')
324 ax.plot(t1, y1, 'r-', label=f'RK4 κ³ μ ({len(t1)} pts)')
325 ax.plot(t2, y2, 'b.', label=f'RKF45 μ μ ({len(t2)} pts)')
326 ax.set_title("μ μμ μ€ν
ν¬κΈ°")
327 ax.legend()
328 ax.grid(True)
329
330 # μμ 4: Lotka-Volterra
331 alpha, beta, gamma, delta = 1.5, 1.0, 3.0, 1.0
332
333 def lotka_volterra(t, state):
334 x, y = state
335 dx = alpha * x - beta * x * y
336 dy = delta * x * y - gamma * y
337 return np.array([dx, dy])
338
339 ax = axes[1, 1]
340 t, sol = rk4_system(lotka_volterra, np.array([10, 5]), (0, 15), 0.01)
341
342 ax.plot(t, sol[:, 0], 'b-', label='νΌμμ')
343 ax.plot(t, sol[:, 1], 'r-', label='ν¬μμ')
344 ax.set_title("Lotka-Volterra λͺ¨λΈ")
345 ax.set_xlabel('μκ°')
346 ax.set_ylabel('κ°μ²΄μ')
347 ax.legend()
348 ax.grid(True)
349
350 plt.tight_layout()
351 plt.savefig('/opt/projects/01_Personal/03_Study/Numerical_Simulation/examples/runge_kutta.png', dpi=150)
352 plt.close()
353 print("κ·Έλν μ μ₯: runge_kutta.png")
354
355
356# =============================================================================
357# ν
μ€νΈ
358# =============================================================================
359def main():
360 print("=" * 60)
361 print("λ£½κ²-μΏ ν λ°©λ² (Runge-Kutta Methods)")
362 print("=" * 60)
363
364 # μμ 1: κΈ°λ³Έ ν
μ€νΈ
365 print("\n[μμ 1] dy/dt = y, y(0) = 1, t β [0, 1]")
366 print("-" * 40)
367
368 f = lambda t, y: y
369 exact = lambda t: np.exp(t)
370
371 t, y = rk4(f, 1, (0, 1), 0.1)
372 print(f"RK4 (h=0.1): y(1) = {y[-1]:.10f}")
373 print(f"μ νκ°: e = {exact(1):.10f}")
374 print(f"μ€μ°¨: = {abs(y[-1] - exact(1)):.2e}")
375
376 # μμ 2: Van der Pol μ§λμ
377 print("\n[μμ 2] Van der Pol μ§λμ")
378 print("-" * 40)
379
380 mu = 1.0
381
382 def van_der_pol(t, state):
383 x, v = state
384 return np.array([v, mu * (1 - x**2) * v - x])
385
386 t, sol = rk4_system(van_der_pol, np.array([2, 0]), (0, 20), 0.01)
387 print(f"μ΄κΈ°: x=2, v=0")
388 print(f"t=20: x={sol[-1, 0]:.4f}, v={sol[-1, 1]:.4f}")
389
390 # μμ 3: Lorenz μμ€ν
(μΉ΄μ€μ€)
391 print("\n[μμ 3] Lorenz μμ€ν
(μΉ΄μ€μ€)")
392 print("-" * 40)
393
394 sigma, rho, beta = 10, 28, 8/3
395
396 def lorenz(t, state):
397 x, y, z = state
398 return np.array([
399 sigma * (y - x),
400 x * (rho - z) - y,
401 x * y - beta * z
402 ])
403
404 t, sol = rk4_system(lorenz, np.array([1, 1, 1]), (0, 50), 0.01)
405 print(f"μ΄κΈ°: (1, 1, 1)")
406 print(f"t=50: ({sol[-1, 0]:.4f}, {sol[-1, 1]:.4f}, {sol[-1, 2]:.4f})")
407 print("(μΉ΄μ€μ€ μμ€ν
- μ΄κΈ°κ° λ―Όκ°μ±)")
408
409 # μμ 4: μ μμ RK45
410 print("\n[μμ 4] μ μμ RK45")
411 print("-" * 40)
412
413 f_test = lambda t, y: -2 * t * y
414 exact_test = lambda t: np.exp(-t**2)
415
416 t, y = rkf45(f_test, 1, (0, 3), tol=1e-8)
417 print(f"dy/dt = -2ty, y(0) = 1")
418 print(f"μ€ν
μ: {len(t)}")
419 print(f"y(3) = {y[-1]:.10f}")
420 print(f"μ ν: {exact_test(3):.10f}")
421 print(f"μ€μ°¨: {abs(y[-1] - exact_test(3)):.2e}")
422
423 # λ°©λ² λΉκ΅
424 compare_methods()
425
426 # μκ°ν
427 try:
428 plot_rk_comparison()
429 except Exception as e:
430 print(f"κ·Έλν μμ± μ€ν¨: {e}")
431
432 print("\n" + "=" * 60)
433 print("λ£½κ²-μΏ ν λ°©λ² μ 리")
434 print("=" * 60)
435 print("""
436 | λ°©λ² | μ°¨μ | ν¨μ νΈμΆ/μ€ν
| νΉμ§ |
437 |----------|------|---------------|-------------------------|
438 | RK2 | 2 | 2 | μμ μ€μΌλ¬μ λμΌ |
439 | RK4 | 4 | 4 | κ°μ₯ λ리 μ¬μ©, μ ν/ν¨μ¨μ |
440 | RK45 | 4(5) | 6 | μ μμ , μ€μ°¨ μΆμ |
441 | RK8 | 8 | 13 | λ§€μ° λμ μ νλ |
442
443 RK4 λΆνλ¬ ν
μ΄λΈ:
444 βββββββ¬ββββββ¬ββββββ¬ββββββ¬ββββββ
445 β 0 β β β β β
446 β 1/2 β 1/2 β β β β
447 β 1/2 β 0 β 1/2 β β β
448 β 1 β 0 β 0 β 1 β β
449 βββββββΌββββββΌββββββΌββββββΌββββββ€
450 β β 1/6 β 1/3 β 1/3 β 1/6 β
451 βββββββ΄ββββββ΄ββββββ΄ββββββ΄ββββββ
452
453 μ€λ¬΄:
454 - scipy.integrate.solve_ivp: λ€μν RK λ°©λ² μ§μ
455 - scipy.integrate.odeint: LSODA (μ μμ λ€λ¨κ³)
456 """)
457
458
459if __name__ == "__main__":
460 main()