import matplotlib.pyplot as plt
from numpy import diff, finfo, double
from ode45 import ode45

def ode45stiff():
	"""tries to solve the logistic differential equation using ode45 with
	adaptive stepsize control and draws the results:
	resulting data points -> logode451.eps
	timestep size         -> logode452.eps
	both plots also contain the 'exact' solution
	"""
	
	# define the differential equation
	fun = lambda t,x: 500*x**2*(1-x)
	# define the time span
	tspan = (0.0,1.0)
	# get total time span
	L = diff(tspan)
	# starting value
	y0 = 0.01
	
	# make the figure
	fig = plt.figure()
	
	# exact solution 
	# set options
	eps = finfo(double).eps
	options = {'reltol':10*eps,'abstol':eps,'stats':'on'}
	tex,yex = ode45(fun,(0,1),y0,**options)
	# plot 'exact' results
	plt.plot(tex,yex,'--',color=[0,0.75,0],linewidth=2,label='$y(t)$')
	plt.hold(True)
	
	# ODE45
	# set options
	options = {'reltol':0.01,'abstol':0.001,'stats':'on'}
	# get the solution using ode45
	[t,y] = ode45(fun,(0,1),y0,**options)
	# plot the solution
	plt.plot(t,y,'r+',label='ode45')
	
	# set label, legend, ....
	plt.xlabel('t',fontsize=14)
	plt.ylabel('y')
	plt.title('ode45 for d_ty = 500 y^2(1-y)')
	plt.show()
	# write to file
	#	plt.savefig('../PICTURES/logode451.eps')
	
	# now plot stepsizes
	fig = plt.figure()
	plt.title('ode45 for d_ty = 500 y^2(1-y)')
	ax1 = fig.add_subplot(111)
	ax1.plot(tex,yex[:,2], 'r--', linewidth=2)
	ax1.set_xlabel('t')
	ax1.set_ylabel('y(t)')
	
	ax2 = ax1.twinx()
	ax2.plot(0.5*(t[:-1]+t[1:]),diff(t), 'r-')
	ax2.set_ylabel('stepsize')
	
	plt.show()
	
	# write to file
	#plt.savefig('../PICTURES/logode452.eps')


if __name__ == '__main__':
	print 'warning: this function takes a long time to complete! (not really worth it)\n'
	ode45stiff()