03_ode_euler.py - Examples

  1"""
  2상미분방정식 - 오일러 방법 (Euler Method)
  3Ordinary Differential Equations - Euler Method
  4
  5초기값 문제 dy/dx = f(x, y), y(x₀) = y₀ 를 수치적으로 푸는 가장 기본적인 방법입니다.
  6"""
  7
  8import numpy as np
  9import matplotlib.pyplot as plt
 10from typing import Callable, Tuple, List
 11
 12
 13# =============================================================================
 14# 1. 전진 오일러 방법 (Forward Euler Method)
 15# =============================================================================
 16def euler_forward(
 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    전진 오일러 방법 (명시적 오일러)
 24
 25    y_{n+1} = y_n + h * f(t_n, y_n)
 26
 27    오차: O(h) - 1차 정확도
 28    안정성: 조건부 안정
 29
 30    Args:
 31        f: dy/dt = f(t, y)
 32        y0: 초기값 y(t0)
 33        t_span: (t0, tf) 시간 구간
 34        h: 시간 간격
 35
 36    Returns:
 37        (t 배열, y 배열)
 38    """
 39    t0, tf = t_span
 40    n_steps = int((tf - t0) / h)
 41
 42    t = np.linspace(t0, tf, n_steps + 1)
 43    y = np.zeros(n_steps + 1)
 44    y[0] = y0
 45
 46    for i in range(n_steps):
 47        y[i + 1] = y[i] + h * f(t[i], y[i])
 48
 49    return t, y
 50
 51
 52# =============================================================================
 53# 2. 후진 오일러 방법 (Backward Euler Method)
 54# =============================================================================
 55def euler_backward(
 56    f: Callable[[float, float], float],
 57    df_dy: Callable[[float, float], float],
 58    y0: float,
 59    t_span: Tuple[float, float],
 60    h: float,
 61    newton_tol: float = 1e-10,
 62    max_iter: int = 50
 63) -> Tuple[np.ndarray, np.ndarray]:
 64    """
 65    후진 오일러 방법 (암시적 오일러)
 66
 67    y_{n+1} = y_n + h * f(t_{n+1}, y_{n+1})
 68
 69    암시적 방정식을 뉴턴법으로 해결
 70    오차: O(h) - 1차 정확도
 71    안정성: 무조건 안정 (stiff 문제에 적합)
 72
 73    Args:
 74        f: dy/dt = f(t, y)
 75        df_dy: ∂f/∂y (야코비안)
 76        y0: 초기값
 77        t_span: 시간 구간
 78        h: 시간 간격
 79    """
 80    t0, tf = t_span
 81    n_steps = int((tf - t0) / h)
 82
 83    t = np.linspace(t0, tf, n_steps + 1)
 84    y = np.zeros(n_steps + 1)
 85    y[0] = y0
 86
 87    for i in range(n_steps):
 88        # 뉴턴법으로 y_{n+1} 구하기
 89        # g(y) = y - y_n - h*f(t_{n+1}, y) = 0
 90        y_new = y[i]  # 초기 추정값 (전진 오일러)
 91        t_new = t[i + 1]
 92
 93        for _ in range(max_iter):
 94            g = y_new - y[i] - h * f(t_new, y_new)
 95            dg = 1 - h * df_dy(t_new, y_new)
 96
 97            if abs(dg) < 1e-15:
 98                break
 99
100            delta = g / dg
101            y_new = y_new - delta
102
103            if abs(delta) < newton_tol:
104                break
105
106        y[i + 1] = y_new
107
108    return t, y
109
110
111# =============================================================================
112# 3. 수정 오일러 방법 (Modified Euler / Heun's Method)
113# =============================================================================
114def euler_modified(
115    f: Callable[[float, float], float],
116    y0: float,
117    t_span: Tuple[float, float],
118    h: float
119) -> Tuple[np.ndarray, np.ndarray]:
120    """
121    수정 오일러 방법 (Heun's Method)
122
123    예측: y*_{n+1} = y_n + h * f(t_n, y_n)
124    수정: y_{n+1} = y_n + h/2 * [f(t_n, y_n) + f(t_{n+1}, y*_{n+1})]
125
126    오차: O(h²) - 2차 정확도
127    RK2의 일종
128    """
129    t0, tf = t_span
130    n_steps = int((tf - t0) / h)
131
132    t = np.linspace(t0, tf, n_steps + 1)
133    y = np.zeros(n_steps + 1)
134    y[0] = y0
135
136    for i in range(n_steps):
137        k1 = f(t[i], y[i])
138        y_pred = y[i] + h * k1
139        k2 = f(t[i + 1], y_pred)
140        y[i + 1] = y[i] + h * (k1 + k2) / 2
141
142    return t, y
143
144
145# =============================================================================
146# 4. 연립 ODE 풀기
147# =============================================================================
148def euler_system(
149    f: Callable[[float, np.ndarray], np.ndarray],
150    y0: np.ndarray,
151    t_span: Tuple[float, float],
152    h: float
153) -> Tuple[np.ndarray, np.ndarray]:
154    """
155    연립 ODE를 전진 오일러로 풀기
156
157    dy/dt = f(t, y), y = [y1, y2, ..., yn]
158
159    Args:
160        f: 벡터 함수 f(t, y) -> dy/dt
161        y0: 초기값 벡터
162        t_span: 시간 구간
163        h: 시간 간격
164
165    Returns:
166        (t 배열, y 배열) - y.shape = (n_steps+1, n_vars)
167    """
168    t0, tf = t_span
169    n_steps = int((tf - t0) / h)
170    n_vars = len(y0)
171
172    t = np.linspace(t0, tf, n_steps + 1)
173    y = np.zeros((n_steps + 1, n_vars))
174    y[0] = y0
175
176    for i in range(n_steps):
177        y[i + 1] = y[i] + h * np.array(f(t[i], y[i]))
178
179    return t, y
180
181
182# =============================================================================
183# 5. 2차 ODE를 1차 연립으로 변환
184# =============================================================================
185def solve_second_order(
186    f: Callable[[float, float, float], float],
187    y0: float,
188    v0: float,
189    t_span: Tuple[float, float],
190    h: float
191) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
192    """
193    2차 ODE: y'' = f(t, y, y')
194    변환: y₁ = y, y₂ = y'
195          y₁' = y₂
196          y₂' = f(t, y₁, y₂)
197
198    Args:
199        f: y'' = f(t, y, y')
200        y0: 초기 위치
201        v0: 초기 속도
202        t_span, h: 시간 구간과 간격
203
204    Returns:
205        (t, y, v) - 위치와 속도
206    """
207    def system(t, state):
208        y, v = state
209        return np.array([v, f(t, y, v)])
210
211    t, solution = euler_system(system, np.array([y0, v0]), t_span, h)
212    return t, solution[:, 0], solution[:, 1]
213
214
215# =============================================================================
216# 오차 분석
217# =============================================================================
218def analyze_euler_error():
219    """오일러 방법의 오차 분석"""
220    # dy/dt = y, y(0) = 1  →  y = e^t
221    f = lambda t, y: y
222    df_dy = lambda t, y: 1
223    y0 = 1
224    t_span = (0, 1)
225    exact = lambda t: np.exp(t)
226
227    print("\n오일러 방법 오차 분석 (dy/dt = y, y(0) = 1)")
228    print("-" * 70)
229    print(f"{'h':>10} | {'전진 오일러':>15} | {'수정 오일러':>15} | {'후진 오일러':>15}")
230    print("-" * 70)
231
232    hs = [0.1, 0.05, 0.025, 0.0125, 0.00625]
233
234    for h in hs:
235        t1, y1 = euler_forward(f, y0, t_span, h)
236        t2, y2 = euler_modified(f, y0, t_span, h)
237        t3, y3 = euler_backward(f, df_dy, y0, t_span, h)
238
239        error1 = abs(y1[-1] - exact(1))
240        error2 = abs(y2[-1] - exact(1))
241        error3 = abs(y3[-1] - exact(1))
242
243        print(f"{h:>10.5f} | {error1:>15.2e} | {error2:>15.2e} | {error3:>15.2e}")
244
245
246# =============================================================================
247# 시각화
248# =============================================================================
249def plot_euler_comparison():
250    """오일러 방법 비교 시각화"""
251    # dy/dt = -2y + sin(t), y(0) = 1
252    f = lambda t, y: -2*y + np.sin(t)
253    df_dy = lambda t, y: -2
254    y0 = 1
255    t_span = (0, 5)
256
257    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
258
259    # 정확해 (scipy 사용)
260    from scipy.integrate import odeint
261    t_exact = np.linspace(0, 5, 500)
262    y_exact = odeint(lambda y, t: f(t, y), y0, t_exact).flatten()
263
264    # 다양한 h 값 비교
265    hs = [0.5, 0.25, 0.1]
266    colors = ['r', 'g', 'b']
267
268    ax = axes[0, 0]
269    ax.plot(t_exact, y_exact, 'k-', linewidth=2, label='정확해')
270    for h, c in zip(hs, colors):
271        t, y = euler_forward(f, y0, t_span, h)
272        ax.plot(t, y, f'{c}o-', markersize=4, label=f'h={h}')
273    ax.set_title('전진 오일러')
274    ax.legend()
275    ax.grid(True)
276
277    ax = axes[0, 1]
278    ax.plot(t_exact, y_exact, 'k-', linewidth=2, label='정확해')
279    for h, c in zip(hs, colors):
280        t, y = euler_modified(f, y0, t_span, h)
281        ax.plot(t, y, f'{c}o-', markersize=4, label=f'h={h}')
282    ax.set_title('수정 오일러')
283    ax.legend()
284    ax.grid(True)
285
286    # 안정성 비교 (stiff 문제)
287    # dy/dt = -50(y - cos(t)), y(0) = 0
288    f_stiff = lambda t, y: -50*(y - np.cos(t))
289    df_stiff = lambda t, y: -50
290    h = 0.05
291    t_span_stiff = (0, 1)
292
293    ax = axes[1, 0]
294    t_ex = np.linspace(0, 1, 500)
295    y_ex = odeint(lambda y, t: f_stiff(t, y), 0, t_ex).flatten()
296    ax.plot(t_ex, y_ex, 'k-', linewidth=2, label='정확해')
297
298    t1, y1 = euler_forward(f_stiff, 0, t_span_stiff, h)
299    t2, y2 = euler_backward(f_stiff, df_stiff, 0, t_span_stiff, h)
300
301    ax.plot(t1, y1, 'r.-', label=f'전진 오일러 (h={h})')
302    ax.plot(t2, y2, 'b.-', label=f'후진 오일러 (h={h})')
303    ax.set_title('Stiff 문제 안정성')
304    ax.legend()
305    ax.grid(True)
306
307    # 조화 진동자
308    # y'' = -y  →  y' = v, v' = -y
309    def harmonic(t, state):
310        y, v = state
311        return np.array([v, -y])
312
313    ax = axes[1, 1]
314    t_ho = np.linspace(0, 20, 1000)
315    y_ho_exact = np.cos(t_ho)  # y(0)=1, y'(0)=0
316
317    t, sol = euler_system(harmonic, np.array([1, 0]), (0, 20), 0.1)
318    ax.plot(t_ho, y_ho_exact, 'k-', linewidth=2, label='정확해')
319    ax.plot(t, sol[:, 0], 'r-', label='오일러 (h=0.1)')
320    ax.set_title('조화 진동자 y\'\' = -y')
321    ax.set_xlabel('t')
322    ax.legend()
323    ax.grid(True)
324
325    plt.tight_layout()
326    plt.savefig('/opt/projects/01_Personal/03_Study/Numerical_Simulation/examples/ode_euler.png', dpi=150)
327    plt.close()
328    print("그래프 저장: ode_euler.png")
329
330
331# =============================================================================
332# 테스트
333# =============================================================================
334def main():
335    print("=" * 60)
336    print("상미분방정식 - 오일러 방법")
337    print("=" * 60)
338
339    # 예제 1: 지수 감쇠 dy/dt = -y
340    print("\n[예제 1] 지수 감쇠: dy/dt = -y, y(0) = 1")
341    print("-" * 40)
342
343    f = lambda t, y: -y
344    exact = lambda t: np.exp(-t)
345
346    t, y_forward = euler_forward(f, 1, (0, 2), 0.2)
347    t, y_modified = euler_modified(f, 1, (0, 2), 0.2)
348
349    print(f"t=2에서:")
350    print(f"  정확값:       {exact(2):.6f}")
351    print(f"  전진 오일러:  {y_forward[-1]:.6f}")
352    print(f"  수정 오일러:  {y_modified[-1]:.6f}")
353
354    # 예제 2: 로지스틱 방정식
355    print("\n[예제 2] 로지스틱 방정식: dy/dt = y(1-y), y(0) = 0.1")
356    print("-" * 40)
357
358    f_logistic = lambda t, y: y * (1 - y)
359    exact_logistic = lambda t: 0.1 * np.exp(t) / (1 + 0.1 * (np.exp(t) - 1))
360
361    t, y = euler_modified(f_logistic, 0.1, (0, 5), 0.1)
362    print(f"t=5에서:")
363    print(f"  정확값:       {exact_logistic(5):.6f}")
364    print(f"  수정 오일러:  {y[-1]:.6f}")
365
366    # 예제 3: 진자 운동 (2차 ODE)
367    print("\n[예제 3] 단진자: θ'' = -sin(θ), θ(0) = π/4, θ'(0) = 0")
368    print("-" * 40)
369
370    f_pendulum = lambda t, theta, omega: -np.sin(theta)
371    t, theta, omega = solve_second_order(f_pendulum, np.pi/4, 0, (0, 10), 0.01)
372
373    print(f"주기적 운동 시뮬레이션 완료")
374    print(f"  초기 각도: {np.degrees(np.pi/4):.1f}°")
375    print(f"  t=10에서 각도: {np.degrees(theta[-1]):.2f}°")
376
377    # 예제 4: Lotka-Volterra (포식자-피식자 모델)
378    print("\n[예제 4] Lotka-Volterra: 포식자-피식자 모델")
379    print("-" * 40)
380
381    alpha, beta, gamma, delta = 1.5, 1.0, 3.0, 1.0
382
383    def lotka_volterra(t, state):
384        x, y = state  # x=피식자, y=포식자
385        dx = alpha * x - beta * x * y
386        dy = delta * x * y - gamma * y
387        return np.array([dx, dy])
388
389    t, sol = euler_system(lotka_volterra, np.array([10, 5]), (0, 15), 0.001)
390    print(f"  초기: 피식자={10}, 포식자={5}")
391    print(f"  t=15: 피식자={sol[-1, 0]:.2f}, 포식자={sol[-1, 1]:.2f}")
392
393    # 오차 분석
394    analyze_euler_error()
395
396    # 시각화
397    try:
398        plot_euler_comparison()
399    except Exception as e:
400        print(f"그래프 생성 실패: {e}")
401
402    print("\n" + "=" * 60)
403    print("오일러 방법 정리")
404    print("=" * 60)
405    print("""
406    | 방법        | 정확도 | 안정성    | 특징                    |
407    |------------|--------|----------|-------------------------|
408    | 전진 오일러 | O(h)   | 조건부   | 가장 단순, 명시적        |
409    | 후진 오일러 | O(h)   | 무조건   | 암시적, stiff에 적합     |
410    | 수정 오일러 | O(h²)  | 조건부   | 2차 정확도, RK2의 일종   |
411
412    한계:
413    - 정확도가 낮음 (더 높은 차수 필요시 RK4 사용)
414    - 에너지 보존 불가 (심플렉틱 적분기 필요)
415
416    실무:
417    - scipy.integrate.odeint: 적응적 다단계 방법
418    - scipy.integrate.solve_ivp: 다양한 방법 선택 가능
419    """)
420
421
422if __name__ == "__main__":
423    main()