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 y, relerror, abserror def rhs(y,t): return -2*y*t ts = N.arange(0,2,0.1) y0 = N.exp(ts[0]**2) true_y = N.exp(-ts**2) ys=odeint( rhs, y0, ts) import pylab for tol,style in zip([1e-6,1e-7],['x','o']): rtol=atol=tol ys=odeint( rhs, y0, ts, atol=atol, rtol=rtol) pylab.subplot(2,1,1) pylab.plot(ts,ys[:,0],style,label='result odeint (rtol=%g, atol=%g)' % (rtol,atol)) pylab.legend() pylab.subplot(2,1,2) pylab.plot(ts,ys[:,0]-true_y,style,label='error (rtol=%g, atol=%g)' % (rtol,atol)) pylab.legend() pylab.xlabel('t') pylab.subplot(2,1,1) pylab.xlabel('t') pylab.plot(ts,true_y,label='exact solution') pylab.legend() pylab.title('Comparing numerical approximations with the exact solution') pylab.savefig('odeintcompareexact.png') pylab.savefig('odeintcompareexact.eps') pylab.show()