PROBLEMS

7.2

For a simple harmonic oscillator $y¨+y = 0$, with initial conditions $y(0) = 1, y˙(0) = 0$, find $y(t)$ from $t = 0$ to 100π. Use an Euler method and a fixed-step fourth-order Runge–Kutta method. For each method check how the average local error, and the error in the final value and slope, depend on the step size.

\begin{equation} \begin{aligned} y(x+h) &= y(x) + h f(x, y(x)) \\ y(x+h) &= y(x) + h \dot y(x) \\ \end{aligned} \end{equation}

Similarly

\begin{equation} \begin{aligned} \ddot y(x+h) &= \dot y(x) + h \ddot y(x) \\ \ddot y(x+h) &= \dot y(x) + h (-y(x)) \end{aligned} \end{equation}

import numpy as np
import time
import matplotlib.pyplot as plt
import math
%matplotlib inline

plt.rcParams['figure.figsize'] = 20, 5

t_start = 0
t_end = 100 * math.pi
h = 0.001

def euler(t_start, t_end, h):
    y_prev = 1
    dy_prev = 0
    
    y = [y_prev]
    steps = [0]
    
    for i in np.arange(t_start, t_end, h):
        y_new = y_prev + h*dy_prev
        dy_new = dy_prev - h*y_prev
        
        y_prev = y_new
        dy_prev = dy_new
        
        y.append(y_new)
        steps.append(i)
        
    return (np.array(y), np.array(steps))

(y_vals, steps) = euler(t_start, t_end, h)

plt.plot(steps, y_vals)
plt.show()