import numpy as np
from numpy import *
from matplotlib import pyplot as plt

# ddy + y = 0

# x1 = y
# x2 = dy

# dx1 = x2 		( = dy)
# dx2 = -x1 	( = ddy)

def f(x1,x2):
	return [x2, -x1]

# Euler:
# y(x + h) = y(x) + h*f(x,y(x))

# steps = [0.001, 0.005, 0.01, 0.02, 0.05]
steps = linspace(0.001,0.1,50)
g_errors = []
l_errors = []
for h in steps:
	[x1, x2] = [1, 0]
	t = arange(0,100*pi,h)
	xx1 = []
	xx2 = []
	x1_exacts = []
	errors = []
	for i in t:
		[dx1, dx2] = f(x1,x2)
		[x1, x2] = [x1 + h*dx1, x2 + h*dx2]
		xx1.append(x1)
		xx2.append(x2)
		x1_exact = cos(i)
		x1_exacts.append(x1_exact)
		error = abs(x1-x1_exact)
		errors.append(error)

	# plt.figure()
	# plt.plot(t,xx1)
	# plt.plot(t,x1_exacts)
	# plt.figure()
	# plt.plot(t,errors)

	print("Step size = %.3f" %(h))
	print("Global error = %.2f" % (error))
	print("Average local error = %.2f" % (mean(errors)))
	g_errors.append(error)
	l_errors.append(mean(errors))

plt.figure()
gerr, = plt.semilogy(steps,g_errors)
lerr, = plt.semilogy(steps,l_errors)
plt.xlabel('Step Size')
plt.ylabel('Error')
plt.legend([gerr, lerr],['Global Error', 'Local Error'],'best')
plt.title('Euler Error vs Step-Size')
plt.show()

