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 absolute local error rtol = 0.0; 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) #test limiting relative local error rtol = 1e-8; atol = 0 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) #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)