import numpy as N
from scipy.integrate import odeint

def testodeint( f, y0, ts, true_y, myrtol, myatol):
    """integrate dy/dt=f(y,t) and compare with exact solution in ture_y.
    Returns the maximum global relative error and
    maximum global absolute error. """
    
    y = odeint( f, y0, ts, rtol=myrtol, atol=myatol )

    abserror = max(abs(true_y-y[:,0]))
    relerror = max(abs((true_y-y[:,0])/true_y))

    return relerror, abserror

def rhs(y,t):
    return -2*y*t

ts = N.arange(0,2,0.01)
y0 = N.exp(ts[0]**2)
true_y = N.exp(-ts**2)

#test limiting one of rel and ans error
rtol = 1e-8; atol = 1e-8
relerr, abserr = testodeint( rhs, y0, ts, true_y, rtol, atol)
print "rtol=%5g, atol=%5g => glob.rel.err=%5g, abs.err=%5g" %(rtol,atol,relerr,abserr)


y = odeint( rhs, y0, ts )
abserr = max(abs(true_y-y[:,0]))
relerr = max(abs((true_y-y[:,0])/true_y))
print "rtol=None, atol=None => glob.rel.err=%5g, abs.err=%5g" %(relerr,abserr)