# Ordinary Different Equation Solving via Runge-Kutta 4

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

%pylab inline


1. (a) Implement RK4.

def _step_rk4(f, t0, y0, dt):
"""
Perform one Fourth-Order Runge-Kutta Step

:param f: function
ODE function with input f(t, y)
:param t0: float
Initial starting time
:param y0: numpy.ndarray()
Initial ODE state (vector)
:param dt: float
Timestep
"""
k1 = f(t0, y0)
k2 = f(t0 + (dt / 2), y0 + (k1 * dt / 2))
k3 = f(t0 + (dt / 2), y0 + (k2 * dt / 2))
k4 = f(t0 + dt, y0 + (k3 * dt))
y = y0 + dt * ((k1 / 6) + (k2 / 3) +
(k3 / 3) + (k4 / 6))
return (t0 + dt,
(y0 + dt * ((k1 / 6) + (k2 / 3) +
(k3 / 3) + (k4 / 6))))

def mrk4(f, t0, y0, dt, n, writecsv=''):
"""
Fixed-Step Fourth-Order Runge-Kutta ODE Solver

:param f: function
ODE function with input f(t, y)
:param t0: float
Initial starting time
:param y0: numpy.ndarray()
Initial ODE state (vector)
:param dt: float
Timestep
:param n: int
Number of iterations (steps) to perform
:param writecsv: bool
:default: False
Write to csv file?
"""
dim = y0.size

# Establish blank solution trajectory
# [[y00, ..., y0n, t0],
#  [y10, ..., y1n, t1],
# ...]
traj = np.zeros((n + 1, dim + 1), dtype=np.float64)

# Set initial position
traj[0, 0:dim] = y0
traj[0, -1]    = t0

# Iterate
for i in range(1, n + 1):
(traj[i, -1],
traj[i, 0:dim]) = _step_rk4(f,
traj[i - 1, -1],
traj[i - 1, 0:dim],
dt)

if writecsv != '':
with open(writecsv, 'w') as f:
csvwriter = csv.writer(f)
[csvwriter.writerow(line) for line in traj]
return traj

k = 1
m = 0.5
garr = [0, 9.8]

fig, axarr = plt.subplots(len(garr), figsize=(10, 10))
for i in range(len(garr)):
g = garr[i]
spring_system = lambda t, y: np.array([y[1],
-(k / m) * y[0] + g])
initial_condition = np.array([1, 0])

rk_points = mrk4(spring_system, 0, initial_condition, 0.1, 100)

axarr[i].plot(rk_points[:, 0], rk_points[:, 1], label=r'$g = {}$'.format(g))
axarr[i].set_title(r'$g = {}$'.format(g))
plt.show()


(b) Compare your results to the solution generated by the forward Euler solution. Which is more accurate?

def forward_euler(f, t0, y0, dt, n):
dim = y0.size

# Establish blank solution trajectory
# [[y00, ..., y0n, t0],
#  [y10, ..., y1n, t1],
# ...]
traj = np.zeros((n + 1, dim + 1), dtype=np.float64)

# Set initial position
traj[0, 0:dim] = y0
traj[0, -1]    = t0

# Iterate
for i in range(1, n + 1):
traj[i, -1] = traj[i - 1, -1] + dt
traj[i, 0:dim] = traj[i - 1, 0:dim] + dt * f(traj[i - 1, -1],
traj[i - 1, 0:dim])
return traj
def reverse_euler(f, t0, y0, dt, n):
dim = y0.size

# Establish blank solution trajectory
# [[y00, ..., y0n, t0],
#  [y10, ..., y1n, t1],
# ...]
traj = np.zeros((n + 1, dim + 1), dtype=np.float64)

# Set initial position
traj[0, 0:dim] = y0
traj[0, -1]    = t0

for i in range(1, n + 1):
traj[i, -1] = traj[i - 1, -1] + dt
y = traj[i - 1, 0:dim]
for j in range(1000):
y = traj[i - 1, 0:dim] + dt * f(traj[i, -1], y)
traj[i, 0:dim] = traj[i - 1, 0:dim] + dt * f(traj[i, -1], y)
return traj
def trapezoid_method(f, t0, y0, dt, n):
dim = y0.size

# Establish blank solution trajectory
# [[y00, ..., y0n, t0],
#  [y10, ..., y1n, t1],
# ...]
traj = np.zeros((n + 1, dim + 1), dtype=np.float64)

# Set initial position
traj[0, 0:dim] = y0
traj[0, -1]    = t0

for i in range(1, n + 1):
traj[i, -1] = traj[i - 1, -1] + dt
traj[i, 0:dim] = (traj[i - 1, 0:dim] +
(dt / 2) *
(f(traj[i - 1, -1], traj[i - 1, 0:dim]) +
f(traj[i - 1, -1] + dt,
traj[i - 1, 0:dim] + dt * f(traj[i - 1, -1],
traj[i - 1, 0:dim]))))
return traj

k = 1
m = 0.5
garr = [0]

for g in garr:
spring_system = lambda t, y: np.array([y[1],
-(k / m) * y[0] + g])
initial_condition = np.array([1, 0])
rk_points = mrk4(spring_system, 0, initial_condition, 0.1, 100)
forward_points = forward_euler(spring_system, 0, initial_condition, 0.1, 100)
reverse_points = reverse_euler(spring_system, 0, initial_condition, 0.1, 100)
trapezoid_points = trapezoid_method(spring_system, 0, initial_condition, 0.1, 100)

fig, axarr = plt.subplots(2, 2, figsize=(10, 10))
axarr[0, 0].plot(rk_points[:, 0], rk_points[:, 1], label='Runge-Kutta')
axarr[0, 1].plot(forward_points[:, 0], forward_points[:, 1], label='Forward Euler')
axarr[1, 0].plot(reverse_points[:, 0], reverse_points[:, 1], label='Reverse Euler')
axarr[1, 1].plot(trapezoid_points[:, 0], trapezoid_points[:, 1], label='Trapezoid Method')
axarr[0, 0].set_title('Runge-Kutta')
axarr[0, 1].set_title('Forward Euler')
axarr[1, 0].set_title('Reverse Euler')
axarr[1, 1].set_title('Trapezoid Method')
plt.suptitle(r'$g = {}$'.format(g))
plt.show()


We can see that Runge-Kutta is far more accurate than forward euler, which spins inward.

(c) Comment on runtime differences

Forward Euler is much faster, as it has less computations to run. This makes much sense.

2. Problem 4(a) and (b) on page 291.

Find the solutions of the IVP given by $y(0) = 0$ and the following first order linear differential equations.

(a) $y^\prime = t + y$

$y(t) = -t + e^t - 1$

(b) $y^\prime = t - y$

$y(t) = t+e^{-t}-1$

4. When and why is it a good idea to use an implicit ODE solver?

One should use an implicit ODE solver when you need to keep the error bounds small, and you have more cycles that you can throw at the program. Because these solvers take more time they are generally implemented less often.