from solvers import *
from pylab import *

# Our system of differential equations
f = (
    lambda t, y: y[1],         # y1' = y2
    lambda t, y: -1.0 * y[0]   # y2' = -y1
)

# Initial conditions
f0 = (1.0, 0.0)

# Step size
steps = [0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.1, 0.05, 0.025, 0.0125]

# Interval
interval = (0.0, 100.0 * pi)

errs_euler = []
errs_rk4 = []

global_errs_euler_value = []
global_errs_rk4_value = []
global_errs_euler_slope = []
global_errs_rk4_slope = []

for h in steps:
    print h
    t1, y1, err_local_euler = solver(euler_step, f, f0, interval, h)
    t2, y2, err_local_rk4 = solver(RK4_step, f, f0, interval, h)
    errs_euler.append(mean(err_local_euler))
    errs_rk4.append(mean(err_local_rk4))

    global_errs_euler_value.append( abs(cos(t1[-1]) - y1[0][-1]) )
    global_errs_rk4_value.append( abs(cos(t2[-1]) - y2[0][-1]) )
    global_errs_euler_slope.append( abs(-sin(t1[-1]) - y1[1][-1]) )
    global_errs_rk4_slope.append( abs(-sin(t2[-1]) - y2[1][-1]) )

"""
subplot(121)
plot(steps, errs_euler)
xlabel("step size")
ylabel("error")
title("Average local error, Euler method")

subplot(122)
plot(steps, errs_rk4)
xlabel("step size")
ylabel("error")
title("Average local error, RK4 method")
"""

"""
subplot(121)
plot(steps, global_errs_euler_value)
xlabel("step size")
ylabel("error")
title("Global value error, Euler method")

subplot(122)
plot(steps, global_errs_rk4_value)
xlabel("step size")
ylabel("error")
title("Global value error, RK4 method")
"""

subplot(121)
plot(steps, global_errs_euler_slope)
xlabel("step size")
ylabel("error")
title("Global slope error, Euler method")

subplot(122)
plot(steps, global_errs_rk4_slope)
xlabel("step size")
ylabel("error")
title("Global slope error, RK4 method")

show()