import scipy.integrate
import matplotlib.pyplot as plt
import numpy
a=1.0
b=2.0
def fun(t):
if t<=-b:
return -a
elif f<b:
return t*a/b
else:
return a
g=lambda t:fun(t)
N=100
time_step=0.1
time_end=10.0
t0=0.0
x0=[[0.5*k,0.5*k] for k in range(-10,10)]
def f(x,t):
return [x[1],-x[0]+g(x[0])]
time_range=[t0..time_end, step=time_step]
plt.figure()
for n in range(10):
sol = scipy.integrate.odeint(f,x0[n],time_range)
x = sol[:,0]
y = sol[:,1]
plt.plot(x,y)
plt.savefig('phase.png')
#corresponding vector field
x0,x1=var('x0 x1')
p=plot_vector_field ((x1,-x0+g(x0)),(x0,-8,8),(x1,-8,8))
--
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
