problem3.py

from pylab import *
from solvers 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
errors = [0.0001, 0.00001, 0.000001, 0.0000001, 0.00000001, 0.000000001]

# Interval
interval = (0.0, 100.0 * pi)

stepsizes = []

for e in errors:
    print e
    t, y, err_local, h_adaptive = solver_adaptive(RK4_step, f, f0, interval,
        e, e*2, 0.3, 750.0)
    stepsizes.append(mean(h_adaptive))

title("Adaptive Step Size")
xlabel("Target error")
ylabel("Step size")
semilogx(errors, stepsizes)
show()