Euler method approximation is too accurate

Question

I'm trying to code Euler's approximation method to solve the harmonic oscillator problem, and then compare it on a plot with the exact solution. However, I'm finding that my approximation seems to be too good (even for larger step sizes), especially since I know the error is supposed to increase as time goes on. Am I doing something wrong in my Euler's method code? It's all written below.

import numpy as np
import matplotlib.pyplot as plt

w = 1.0 #frequency of oscillation
x_0 = 1.0 #initial position
h = 0.6 #step size
t_0 = 0 #starting time
n = 20
t = np.arange(t_0, n, h) #t = t_0 + nh (h has units of seconds)
y = np.array([[0.0, 0.0]]*len(t))

def x(t, x_0, w): #exact solution
    return x_0*np.cos(w*t)
def x_prime(t, x_0, w): #first derivative of solution (which is also v)
    return -x_0*w*np.sin(w*t)

y[0] = [x_0, 0.0] #starting velocity is 0
for i in range(1, len(t)):
    # Euler's method
    y[i] = y[i-1] + np.array([x_prime(t[i-1], x_0, w), -w**2 * x(t[i-1], x_0, w)]) * h

xlist = [item[0] for item in y] #list of x values at each t value
vlist = [item[1] for item in y] #list of velocity values at each t value

plt.plot(t, xlist, label='Euler approximation')
plt.plot(t, x(t, x_0, w), label='Exact solution')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend(loc='upper right')

And my plot for this code outputs as:

Answer 1

Your problem is that you compute the slopes for the Euler method from the exact solution. You should use the previous approximation instead,

# Euler's method
y[i] = y[i-1] + np.array([y[i-1,1], -w**2 * y[i-1,0]]) * h

or use a function

def f(t,y): x, dotx = y; return np.array([ dotx, -w**2*x])

...

for i in range(1, len(t)):
    # Euler's method
    y[i] = y[i-1] + f(t[i-1],y[i-1]) * h

Then you will get the expected rapidly diverging behavior, here for h=0.1:

---

Btw., it might be simpler to write

xlist, vlist = y.T