Runge Kutta Method 실습 (1일차)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

h = 0.05                                          # step size
t = np.arange(0, 100+h, h)                        # time
u = np.zeros_like(t)                              # position
u[0] = -0.5                                       # initial condition

def F_ut(u, t):
    return 1 * np.exp(-u) - 0.1 * t               # velocity is dependent on time t and position x

# Runge-Kutta method
for i in range(len(t)-1):                         # calculation loop
    k1 = F_ut(u[i], t[i])
    k2 = F_ut(u[i] + 0.5 * h * k1, t[i] + 0.5 * h)
    k3 = F_ut(u[i] + 0.5 * h * k2, t[i] + 0.5 * h)
    k4 = F_ut(u[i] + k3 * h, t[i] + h)
    u[i+1] = u[i] + (1/6) * (k1 + 2*k2 + 2*k3 + k4) * h

plt.plot(t, u, 'b-', label='Runge-Kutta')

# Simple forward
u_simple = np.zeros_like(t)
u_simple[0] = -0.5                                 # initial condition
for i in range(len(t)-1):                          # calculation loop
    u_simple[i+1] = u_simple[i] + h * F_ut(u_simple[i], t[i])

plt.plot(t, u_simple, 'r-.', label='Simple Forward')

# Validate using a decent ODE integrator
def F_ut2(u, t):
    return 1 * np.exp(-u) - 0.1 * t               # time must be specified first
t_span = [0, 100]
tsol = np.linspace(0, 100, len(t))
usol = odeint(F_ut2, u[0], tsol)

plt.plot(tsol, usol, 'g--', label='ODE Integrator')

plt.xlabel('Time')
plt.ylabel('Position')
plt.legend()
plt.show()