ODE ์ˆ˜์น˜ํ•ด๋ฒ•

Previous ์ƒ๋ฏธ๋ถ„๋ฐฉ์ •์‹ ๊ธฐ์ดˆ
๐Ÿ“‹ Topic list
Next ODE ๊ณ ๊ธ‰

ODE ์ˆ˜์น˜ํ•ด๋ฒ•

๊ฐœ์š”

์ƒ๋ฏธ๋ถ„๋ฐฉ์ •์‹์˜ ์ˆ˜์น˜ํ•ด๋ฒ•์€ ํ•ด์„์  ํ•ด๋ฅผ ๊ตฌํ•  ์ˆ˜ ์—†๊ฑฐ๋‚˜ ์–ด๋ ค์šด ๊ฒฝ์šฐ์— ๊ทผ์‚ฌํ•ด๋ฅผ ๊ณ„์‚ฐํ•ฉ๋‹ˆ๋‹ค. ์˜ค์ผ๋Ÿฌ ๋ฐฉ๋ฒ•๋ถ€ํ„ฐ 4์ฐจ ๋ฃฝ๊ฒŒ-์ฟ ํƒ€๊นŒ์ง€ ์ฃผ์š” ๋ฐฉ๋ฒ•๋“ค์„ ํ•™์Šตํ•ฉ๋‹ˆ๋‹ค.

1. ์˜ค์ผ๋Ÿฌ ๋ฐฉ๋ฒ• (Euler Method)

1.1 ์ „์ง„ ์˜ค์ผ๋Ÿฌ (Forward Euler)

import numpy as np
import matplotlib.pyplot as plt

def forward_euler(f, y0, t_span, n_steps):
    """
    ์ „์ง„ ์˜ค์ผ๋Ÿฌ ๋ฐฉ๋ฒ•

    dy/dt = f(t, y)
    y_{n+1} = y_n + h * f(t_n, y_n)
    """
    t = np.linspace(t_span[0], t_span[1], n_steps + 1)
    h = t[1] - t[0]

    y = np.zeros(n_steps + 1)
    y[0] = y0

    for i in range(n_steps):
        y[i+1] = y[i] + h * f(t[i], y[i])

    return t, y

# ํ…Œ์ŠคํŠธ: dy/dt = -y, y(0) = 1
# ํ•ด์„์  ํ•ด: y = e^(-t)
f = lambda t, y: -y
y0 = 1.0
t_span = (0, 5)

# ๋‹ค์–‘ํ•œ ์Šคํ… ์ˆ˜๋กœ ๋น„๊ต
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

t_exact = np.linspace(0, 5, 200)
y_exact = np.exp(-t_exact)

for n in [10, 20, 50, 100]:
    t, y = forward_euler(f, y0, t_span, n)
    axes[0].plot(t, y, 'o-', markersize=3, label=f'n={n}')

axes[0].plot(t_exact, y_exact, 'k-', linewidth=2, label='์ •ํ™•ํ•ด')
axes[0].set_xlabel('t')
axes[0].set_ylabel('y')
axes[0].set_title('์ „์ง„ ์˜ค์ผ๋Ÿฌ ๋ฐฉ๋ฒ•')
axes[0].legend()
axes[0].grid(True)

# ์˜ค์ฐจ ๋ถ„์„
errors = []
n_values = [10, 20, 50, 100, 200, 500]
for n in n_values:
    t, y = forward_euler(f, y0, t_span, n)
    error = np.abs(y[-1] - np.exp(-5))
    errors.append(error)

axes[1].loglog(n_values, errors, 'bo-', label='์‹ค์ œ ์˜ค์ฐจ')
axes[1].loglog(n_values, [5/n for n in n_values], 'r--', label='O(h)')
axes[1].set_xlabel('์Šคํ… ์ˆ˜ n')
axes[1].set_ylabel('์˜ค์ฐจ')
axes[1].set_title('์˜ค์ฐจ vs ์Šคํ… ์ˆ˜ (O(h) ์ˆ˜๋ ด)')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.show()

1.2 ํ›„์ง„ ์˜ค์ผ๋Ÿฌ (Backward Euler)

def backward_euler(f, df_dy, y0, t_span, n_steps, tol=1e-10):
    """
    ํ›„์ง„ ์˜ค์ผ๋Ÿฌ ๋ฐฉ๋ฒ• (์•”์‹œ์ )

    y_{n+1} = y_n + h * f(t_{n+1}, y_{n+1})

    ๋‰ดํ„ด-๋žฉ์Šจ์œผ๋กœ ์•”์‹œ์  ๋ฐฉ์ •์‹ ํ’€์ด
    """
    t = np.linspace(t_span[0], t_span[1], n_steps + 1)
    h = t[1] - t[0]

    y = np.zeros(n_steps + 1)
    y[0] = y0

    for i in range(n_steps):
        # ๋‰ดํ„ด-๋žฉ์Šจ: g(y_{n+1}) = y_{n+1} - y_n - h*f(t_{n+1}, y_{n+1}) = 0
        y_guess = y[i]  # ์ดˆ๊ธฐ ์ถ”์ธก

        for _ in range(100):  # ์ตœ๋Œ€ ๋ฐ˜๋ณต ํšŸ์ˆ˜
            g = y_guess - y[i] - h * f(t[i+1], y_guess)
            dg = 1 - h * df_dy(t[i+1], y_guess)

            y_new = y_guess - g / dg

            if abs(y_new - y_guess) < tol:
                break
            y_guess = y_new

        y[i+1] = y_new

    return t, y

# ํ…Œ์ŠคํŠธ: dy/dt = -y
f = lambda t, y: -y
df_dy = lambda t, y: -1  # โˆ‚f/โˆ‚y

t_fw, y_fw = forward_euler(f, 1.0, (0, 5), 20)
t_bw, y_bw = backward_euler(f, df_dy, 1.0, (0, 5), 20)
t_exact = np.linspace(0, 5, 200)
y_exact = np.exp(-t_exact)

plt.figure(figsize=(10, 5))
plt.plot(t_fw, y_fw, 'bo-', label='์ „์ง„ ์˜ค์ผ๋Ÿฌ')
plt.plot(t_bw, y_bw, 'rs-', label='ํ›„์ง„ ์˜ค์ผ๋Ÿฌ')
plt.plot(t_exact, y_exact, 'k-', linewidth=2, label='์ •ํ™•ํ•ด')
plt.xlabel('t')
plt.ylabel('y')
plt.title('์ „์ง„ vs ํ›„์ง„ ์˜ค์ผ๋Ÿฌ')
plt.legend()
plt.grid(True)
plt.show()

2. ๋ฃฝ๊ฒŒ-์ฟ ํƒ€ ๋ฐฉ๋ฒ• (Runge-Kutta Methods)

2.1 RK2 (์ค‘์ ๋ฒ•, Heun)

def rk2_midpoint(f, y0, t_span, n_steps):
    """RK2 ์ค‘์ ๋ฒ•"""
    t = np.linspace(t_span[0], t_span[1], n_steps + 1)
    h = t[1] - t[0]

    y = np.zeros(n_steps + 1)
    y[0] = y0

    for i in range(n_steps):
        k1 = f(t[i], y[i])
        k2 = f(t[i] + h/2, y[i] + h/2 * k1)
        y[i+1] = y[i] + h * k2

    return t, y

def rk2_heun(f, y0, t_span, n_steps):
    """RK2 Heun ๋ฐฉ๋ฒ• (์ˆ˜์ • ์˜ค์ผ๋Ÿฌ)"""
    t = np.linspace(t_span[0], t_span[1], n_steps + 1)
    h = t[1] - t[0]

    y = np.zeros(n_steps + 1)
    y[0] = y0

    for i in range(n_steps):
        k1 = f(t[i], y[i])
        k2 = f(t[i] + h, y[i] + h * k1)
        y[i+1] = y[i] + h/2 * (k1 + k2)

    return t, y

# ๋น„๊ต
f = lambda t, y: -y
t_span = (0, 5)
n = 20

t_euler, y_euler = forward_euler(f, 1.0, t_span, n)
t_mid, y_mid = rk2_midpoint(f, 1.0, t_span, n)
t_heun, y_heun = rk2_heun(f, 1.0, t_span, n)
t_exact = np.linspace(0, 5, 200)
y_exact = np.exp(-t_exact)

plt.figure(figsize=(10, 5))
plt.plot(t_euler, y_euler, 'o-', label='Euler', alpha=0.7)
plt.plot(t_mid, y_mid, 's-', label='RK2 ์ค‘์ ๋ฒ•', alpha=0.7)
plt.plot(t_heun, y_heun, '^-', label='RK2 Heun', alpha=0.7)
plt.plot(t_exact, y_exact, 'k-', linewidth=2, label='์ •ํ™•ํ•ด')
plt.xlabel('t')
plt.ylabel('y')
plt.title('์˜ค์ผ๋Ÿฌ vs RK2 ๋น„๊ต')
plt.legend()
plt.grid(True)
plt.show()

2.2 RK4 (๊ณ ์ „์  4์ฐจ)

def rk4(f, y0, t_span, n_steps):
    """
    ๊ณ ์ „์  4์ฐจ ๋ฃฝ๊ฒŒ-์ฟ ํƒ€ ๋ฐฉ๋ฒ•

    k1 = f(t_n, y_n)
    k2 = f(t_n + h/2, y_n + h*k1/2)
    k3 = f(t_n + h/2, y_n + h*k2/2)
    k4 = f(t_n + h, y_n + h*k3)

    y_{n+1} = y_n + h/6 * (k1 + 2*k2 + 2*k3 + k4)
    """
    t = np.linspace(t_span[0], t_span[1], n_steps + 1)
    h = t[1] - t[0]

    y = np.zeros(n_steps + 1)
    y[0] = y0

    for i in range(n_steps):
        k1 = f(t[i], y[i])
        k2 = f(t[i] + h/2, y[i] + h/2 * k1)
        k3 = f(t[i] + h/2, y[i] + h/2 * k2)
        k4 = f(t[i] + h, y[i] + h * k3)

        y[i+1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4)

    return t, y

# ์ •ํ™•๋„ ๋น„๊ต
f = lambda t, y: -y
t_span = (0, 5)

methods = [
    ('Euler', forward_euler),
    ('RK2', rk2_heun),
    ('RK4', rk4)
]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# ํ•ด ๋น„๊ต
n = 10
for name, method in methods:
    t, y = method(f, 1.0, t_span, n)
    axes[0].plot(t, y, 'o-', label=name)

t_exact = np.linspace(0, 5, 200)
axes[0].plot(t_exact, np.exp(-t_exact), 'k-', linewidth=2, label='์ •ํ™•ํ•ด')
axes[0].set_xlabel('t')
axes[0].set_ylabel('y')
axes[0].set_title(f'n={n} ์Šคํ…')
axes[0].legend()
axes[0].grid(True)

# ์ˆ˜๋ ด ์ฐจ์ˆ˜
n_values = [10, 20, 50, 100, 200]
for name, method in methods:
    errors = []
    for n in n_values:
        t, y = method(f, 1.0, t_span, n)
        error = np.abs(y[-1] - np.exp(-5))
        errors.append(error)
    axes[1].loglog(n_values, errors, 'o-', label=name)

# ์ด๋ก ์  ์ˆ˜๋ ด ์ฐจ์ˆ˜
axes[1].loglog(n_values, [1/n for n in n_values], 'k--', alpha=0.5, label='O(h)')
axes[1].loglog(n_values, [1/n**2 for n in n_values], 'k:', alpha=0.5, label='O(hยฒ)')
axes[1].loglog(n_values, [1/n**4 for n in n_values], 'k-.', alpha=0.5, label='O(hโด)')

axes[1].set_xlabel('์Šคํ… ์ˆ˜ n')
axes[1].set_ylabel('์˜ค์ฐจ')
axes[1].set_title('์ˆ˜๋ ด ์ฐจ์ˆ˜ ๋น„๊ต')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.show()

3. ์‹œ์Šคํ…œ์— ๋Œ€ํ•œ ์ ์šฉ

3.1 ๋ฒกํ„ฐํ™”๋œ RK4

def rk4_system(f, y0, t_span, n_steps):
    """
    ์—ฐ๋ฆฝ ODE ์‹œ์Šคํ…œ์„ ์œ„ํ•œ RK4

    y: ๋ฒกํ„ฐ
    f: ๋ฒกํ„ฐ ํ•จ์ˆ˜
    """
    t = np.linspace(t_span[0], t_span[1], n_steps + 1)
    h = t[1] - t[0]

    y = np.zeros((n_steps + 1, len(y0)))
    y[0] = y0

    for i in range(n_steps):
        k1 = np.array(f(t[i], y[i]))
        k2 = np.array(f(t[i] + h/2, y[i] + h/2 * k1))
        k3 = np.array(f(t[i] + h/2, y[i] + h/2 * k2))
        k4 = np.array(f(t[i] + h, y[i] + h * k3))

        y[i+1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4)

    return t, y

# ๋‹จ์ˆœ ์กฐํ™” ์ง„๋™์ž: x'' + ฯ‰ยฒx = 0
# ์‹œ์Šคํ…œ: yโ‚' = yโ‚‚, yโ‚‚' = -ฯ‰ยฒyโ‚
omega = 2.0

def harmonic_oscillator(t, y):
    return [y[1], -omega**2 * y[0]]

t_span = (0, 10)
y0 = [1.0, 0.0]  # x(0) = 1, x'(0) = 0

t, y = rk4_system(harmonic_oscillator, y0, t_span, 200)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# ์‹œ๊ฐ„ ์‘๋‹ต
axes[0].plot(t, y[:, 0], label='x(t)')
axes[0].plot(t, y[:, 1], label="x'(t)")
axes[0].plot(t, np.cos(omega * t), 'k--', alpha=0.5, label='์ •ํ™•ํ•ด')
axes[0].set_xlabel('t')
axes[0].set_ylabel('๊ฐ’')
axes[0].set_title('์กฐํ™” ์ง„๋™์ž')
axes[0].legend()
axes[0].grid(True)

# ์œ„์ƒ ํ‰๋ฉด
axes[1].plot(y[:, 0], y[:, 1])
axes[1].set_xlabel('x')
axes[1].set_ylabel("x'")
axes[1].set_title('์œ„์ƒ ๊ณต๊ฐ„')
axes[1].set_aspect('equal')
axes[1].grid(True)

plt.tight_layout()
plt.show()

3.2 ๊ฐ์‡  ์ง„๋™์ž

def damped_oscillator_example():
    """๊ฐ์‡  ์ง„๋™์ž: x'' + 2ฮณx' + ฯ‰โ‚€ยฒx = 0"""
    omega0 = 2.0
    gamma = 0.3  # ๊ฐ์‡  ๊ณ„์ˆ˜

    def damped_osc(t, y):
        return [y[1], -2*gamma*y[1] - omega0**2 * y[0]]

    t_span = (0, 15)
    y0 = [1.0, 0.0]

    t, y = rk4_system(damped_osc, y0, t_span, 300)

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # ์‹œ๊ฐ„ ์‘๋‹ต
    omega_d = np.sqrt(omega0**2 - gamma**2)  # ๊ฐ์‡  ์ง„๋™์ˆ˜
    envelope = np.exp(-gamma * t)

    axes[0].plot(t, y[:, 0], label='x(t)')
    axes[0].plot(t, envelope, 'r--', label='๊ฐ์‡  envelope')
    axes[0].plot(t, -envelope, 'r--')
    axes[0].set_xlabel('t')
    axes[0].set_ylabel('x')
    axes[0].set_title(f'๊ฐ์‡  ์ง„๋™์ž (ฮณ={gamma}, ฯ‰โ‚€={omega0})')
    axes[0].legend()
    axes[0].grid(True)

    # ์œ„์ƒ ๊ณต๊ฐ„ (๋‚˜์„ ํ˜•)
    axes[1].plot(y[:, 0], y[:, 1])
    axes[1].scatter([0], [0], color='red', s=100, zorder=5, label='์•ˆ์ • ํ‰ํ˜•')
    axes[1].set_xlabel('x')
    axes[1].set_ylabel("x'")
    axes[1].set_title('์œ„์ƒ ๊ณต๊ฐ„')
    axes[1].legend()
    axes[1].grid(True)

    plt.tight_layout()
    plt.show()

damped_oscillator_example()

4. ์ ์‘ํ˜• ์Šคํ… ํฌ๊ธฐ

4.1 ์˜ค์ฐจ ์ถ”์ •

def rk45_step(f, t, y, h):
    """
    RK4-5 (Dormand-Prince) ํ•œ ์Šคํ…

    4์ฐจ์™€ 5์ฐจ ๊ทผ์‚ฌ๋ฅผ ๋™์‹œ์— ๊ณ„์‚ฐํ•˜์—ฌ ์˜ค์ฐจ ์ถ”์ •
    """
    # Dormand-Prince ๊ณ„์ˆ˜
    c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1]
    a = [
        [],
        [1/5],
        [3/40, 9/40],
        [44/45, -56/15, 32/9],
        [19372/6561, -25360/2187, 64448/6561, -212/729],
        [9017/3168, -355/33, 46732/5247, 49/176, -5103/18656],
        [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]
    ]
    b5 = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]
    b4 = [5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]

    # k ๊ณ„์‚ฐ
    k = [np.array(f(t, y))]
    for i in range(1, 7):
        yi = y.copy()
        for j in range(i):
            yi = yi + h * a[i][j] * k[j]
        k.append(np.array(f(t + c[i]*h, yi)))

    # 5์ฐจ ๊ทผ์‚ฌ
    y5 = y.copy()
    for i in range(7):
        y5 = y5 + h * b5[i] * k[i]

    # 4์ฐจ ๊ทผ์‚ฌ (์˜ค์ฐจ ์ถ”์ •์šฉ)
    y4 = y.copy()
    for i in range(7):
        y4 = y4 + h * b4[i] * k[i]

    # ์˜ค์ฐจ ์ถ”์ •
    error = np.max(np.abs(y5 - y4))

    return y5, error

def adaptive_rk45(f, y0, t_span, tol=1e-6, h_init=0.1):
    """์ ์‘ํ˜• RK4-5 ์†”๋ฒ„"""
    t = [t_span[0]]
    y = [np.array(y0)]
    h = h_init

    while t[-1] < t_span[1]:
        # ์Šคํ… ํฌ๊ธฐ๊ฐ€ ๋ฒ”์œ„๋ฅผ ๋ฒ—์–ด๋‚˜์ง€ ์•Š๋„๋ก
        if t[-1] + h > t_span[1]:
            h = t_span[1] - t[-1]

        # RK4-5 ์Šคํ…
        y_new, error = rk45_step(f, t[-1], y[-1], h)

        # ์˜ค์ฐจ์— ๋”ฐ๋ฅธ ์Šคํ… ์กฐ์ ˆ
        if error < tol:
            t.append(t[-1] + h)
            y.append(y_new)

            # ์Šคํ… ํฌ๊ธฐ ์ฆ๊ฐ€ (์•ˆ์ „ ๊ณ„์ˆ˜ 0.9)
            if error > 0:
                h = min(h * 0.9 * (tol / error) ** 0.2, 2 * h)
        else:
            # ์Šคํ… ํฌ๊ธฐ ๊ฐ์†Œ
            h = max(h * 0.9 * (tol / error) ** 0.25, h / 4)

    return np.array(t), np.array(y)

# ํ…Œ์ŠคํŠธ
def stiff_like(t, y):
    """๊ธ‰๊ฒฉํ•œ ๋ณ€ํ™”๊ฐ€ ์žˆ๋Š” ํ•จ์ˆ˜"""
    return [-10 * y[0] if t < 1 else -0.1 * y[0]]

t_adapt, y_adapt = adaptive_rk45(stiff_like, [1.0], (0, 5), tol=1e-6)

plt.figure(figsize=(12, 4))

plt.subplot(1, 2, 1)
plt.plot(t_adapt, y_adapt[:, 0], 'o-', markersize=3)
plt.xlabel('t')
plt.ylabel('y')
plt.title('์ ์‘ํ˜• RK4-5')
plt.grid(True)

plt.subplot(1, 2, 2)
dt = np.diff(t_adapt)
plt.plot(t_adapt[:-1], dt, 'o-')
plt.xlabel('t')
plt.ylabel('์Šคํ… ํฌ๊ธฐ h')
plt.title('์Šคํ… ํฌ๊ธฐ ๋ณ€ํ™”')
plt.grid(True)

plt.tight_layout()
plt.show()

5. SciPy๋ฅผ ์ด์šฉํ•œ ODE ํ’€์ด

5.1 solve_ivp

from scipy.integrate import solve_ivp

# ๋กœ๋ Œ์ธ  ์‹œ์Šคํ…œ
def lorenz(t, state, sigma=10, rho=28, beta=8/3):
    x, y, z = state
    return [
        sigma * (y - x),
        x * (rho - z) - y,
        x * y - beta * z
    ]

t_span = (0, 50)
y0 = [1.0, 1.0, 1.0]

# ๋‹ค์–‘ํ•œ ์†”๋ฒ„ ๋น„๊ต
solvers = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF']

fig = plt.figure(figsize=(15, 10))

for idx, method in enumerate(solvers, 1):
    sol = solve_ivp(lorenz, t_span, y0, method=method,
                    dense_output=True, max_step=0.01)

    ax = fig.add_subplot(2, 3, idx, projection='3d')
    ax.plot(sol.y[0], sol.y[1], sol.y[2], linewidth=0.5)
    ax.set_title(f'{method} (n={len(sol.t)})')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')

plt.tight_layout()
plt.show()

5.2 ์ด๋ฒคํŠธ ๊ฒ€์ถœ

def projectile_motion():
    """ํฌ๋ฌผ์„  ์šด๋™ - ์ง€๋ฉด ์ถฉ๋Œ ๊ฒ€์ถœ"""

    g = 9.8

    def projectile(t, state):
        x, y, vx, vy = state
        return [vx, vy, 0, -g]

    def hit_ground(t, state):
        return state[1]  # y = 0์ผ ๋•Œ ์ด๋ฒคํŠธ

    hit_ground.terminal = True  # ์ด๋ฒคํŠธ ๋ฐœ์ƒ ์‹œ ์ค‘๋‹จ
    hit_ground.direction = -1   # ๊ฐ์†Œํ•  ๋•Œ๋งŒ (๋‚™ํ•˜ ์ค‘)

    # ์ดˆ๊ธฐ ์กฐ๊ฑด: 45๋„ ๊ฐ๋„, ์†๋„ 20 m/s
    v0 = 20
    angle = 45 * np.pi / 180
    y0 = [0, 0, v0 * np.cos(angle), v0 * np.sin(angle)]

    sol = solve_ivp(projectile, (0, 10), y0, events=hit_ground,
                    dense_output=True, max_step=0.01)

    print(f"๋น„ํ–‰ ์‹œ๊ฐ„: {sol.t_events[0][0]:.4f} s")
    print(f"๋น„ํ–‰ ๊ฑฐ๋ฆฌ: {sol.y_events[0][0][0]:.4f} m")

    plt.figure(figsize=(10, 5))
    plt.plot(sol.y[0], sol.y[1])
    plt.scatter(sol.y_events[0][0][0], 0, color='red', s=100, label='์ฐฉ์ง€์ ')
    plt.xlabel('x (m)')
    plt.ylabel('y (m)')
    plt.title('ํฌ๋ฌผ์„  ์šด๋™')
    plt.legend()
    plt.grid(True)
    plt.axis('equal')
    plt.show()

projectile_motion()

6. ์•ˆ์ •์„ฑ ๋ถ„์„

6.1 ์•ˆ์ • ์˜์—ญ

def stability_regions():
    """์ˆ˜์น˜ ๋ฐฉ๋ฒ•์˜ ์•ˆ์ • ์˜์—ญ"""

    # ํ…Œ์ŠคํŠธ ๋ฐฉ์ •์‹: dy/dt = ฮปy
    # ์ •ํ™•ํ•ด: y(nh) = e^(ฮปnh)
    # ์ˆ˜์น˜ํ•ด: y_n = R(ฮปh)^n * y_0

    # ๋ณต์†Œ ํ‰๋ฉด์—์„œ |R(z)| โ‰ค 1์ธ ์˜์—ญ

    x = np.linspace(-4, 2, 200)
    y = np.linspace(-3, 3, 200)
    X, Y = np.meshgrid(x, y)
    Z = X + 1j * Y  # z = ฮปh

    fig, axes = plt.subplots(2, 2, figsize=(12, 10))

    # ์ „์ง„ ์˜ค์ผ๋Ÿฌ: R(z) = 1 + z
    R_euler = np.abs(1 + Z)
    axes[0, 0].contour(X, Y, R_euler, levels=[1], colors='blue')
    axes[0, 0].contourf(X, Y, R_euler, levels=[0, 1], alpha=0.3)
    axes[0, 0].set_title('์ „์ง„ ์˜ค์ผ๋Ÿฌ: |1 + z| โ‰ค 1')
    axes[0, 0].set_xlabel('Re(ฮปh)')
    axes[0, 0].set_ylabel('Im(ฮปh)')
    axes[0, 0].grid(True)
    axes[0, 0].set_aspect('equal')

    # ํ›„์ง„ ์˜ค์ผ๋Ÿฌ: R(z) = 1/(1-z)
    R_backward = np.abs(1 / (1 - Z))
    axes[0, 1].contour(X, Y, R_backward, levels=[1], colors='blue')
    axes[0, 1].contourf(X, Y, R_backward, levels=[0, 1], alpha=0.3)
    axes[0, 1].set_title('ํ›„์ง„ ์˜ค์ผ๋Ÿฌ: |1/(1-z)| โ‰ค 1')
    axes[0, 1].set_xlabel('Re(ฮปh)')
    axes[0, 1].set_ylabel('Im(ฮปh)')
    axes[0, 1].grid(True)
    axes[0, 1].set_aspect('equal')

    # RK2: R(z) = 1 + z + zยฒ/2
    R_rk2 = np.abs(1 + Z + Z**2/2)
    axes[1, 0].contour(X, Y, R_rk2, levels=[1], colors='blue')
    axes[1, 0].contourf(X, Y, R_rk2, levels=[0, 1], alpha=0.3)
    axes[1, 0].set_title('RK2: |1 + z + zยฒ/2| โ‰ค 1')
    axes[1, 0].set_xlabel('Re(ฮปh)')
    axes[1, 0].set_ylabel('Im(ฮปh)')
    axes[1, 0].grid(True)
    axes[1, 0].set_aspect('equal')

    # RK4: R(z) = 1 + z + zยฒ/2 + zยณ/6 + zโด/24
    R_rk4 = np.abs(1 + Z + Z**2/2 + Z**3/6 + Z**4/24)
    axes[1, 1].contour(X, Y, R_rk4, levels=[1], colors='blue')
    axes[1, 1].contourf(X, Y, R_rk4, levels=[0, 1], alpha=0.3)
    axes[1, 1].set_title('RK4')
    axes[1, 1].set_xlabel('Re(ฮปh)')
    axes[1, 1].set_ylabel('Im(ฮปh)')
    axes[1, 1].grid(True)
    axes[1, 1].set_aspect('equal')

    plt.tight_layout()
    plt.show()

stability_regions()

์—ฐ์Šต ๋ฌธ์ œ

๋ฌธ์ œ 1

RK4๋กœ Van der Pol ์ง„๋™์ž๋ฅผ ํ’€์–ด๋ณด์„ธ์š”: x'' - ฮผ(1 - xยฒ)x' + x = 0, ฮผ = 2

def exercise_1():
    mu = 2.0

    def van_der_pol(t, y):
        return [y[1], mu * (1 - y[0]**2) * y[1] - y[0]]

    sol = solve_ivp(van_der_pol, (0, 30), [0.1, 0],
                    method='RK45', max_step=0.01)

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    axes[0].plot(sol.t, sol.y[0])
    axes[0].set_xlabel('t')
    axes[0].set_ylabel('x')
    axes[0].set_title('Van der Pol ์ง„๋™์ž')

    axes[1].plot(sol.y[0], sol.y[1])
    axes[1].set_xlabel('x')
    axes[1].set_ylabel("x'")
    axes[1].set_title('์œ„์ƒ ๊ณต๊ฐ„ (๋ฆฌ๋ฐ‹ ์‚ฌ์ดํด)')

    plt.tight_layout()
    plt.show()

exercise_1()

์š”์•ฝ

๋ฐฉ๋ฒ• ์ฐจ์ˆ˜ ํŠน์ง•
์ „์ง„ ์˜ค์ผ๋Ÿฌ O(h) ๊ฐ„๋‹จ, ์ œํ•œ์  ์•ˆ์ •์„ฑ
ํ›„์ง„ ์˜ค์ผ๋Ÿฌ O(h) ์•”์‹œ์ , A-์•ˆ์ •
RK2 O(hยฒ) ์ค‘์ ๋ฒ•, Heun
RK4 O(hโด) ๊ฐ€์žฅ ๋„๋ฆฌ ์‚ฌ์šฉ
RK4-5 O(hโต) ์ ์‘ํ˜• ์Šคํ…
SciPy ์†”๋ฒ„ ์œ ํ˜• ์šฉ๋„
RK45 ๋ช…์‹œ์  ์ผ๋ฐ˜ ๋ฌธ์ œ (๊ธฐ๋ณธ)
DOP853 ๋ช…์‹œ์  ๊ณ ์ •๋ฐ€๋„
Radau ์•”์‹œ์  ๊ฐ•์„ฑ ๋ฌธ์ œ
BDF ์•”์‹œ์  ๊ฐ•์„ฑ ๋ฌธ์ œ
Previous ์ƒ๋ฏธ๋ถ„๋ฐฉ์ •์‹ ๊ธฐ์ดˆ
๐Ÿ“‹ Topic list
Next ODE ๊ณ ๊ธ‰
← → to navigate between lessons