CoCalc Public Filescloud-examples / sage / 2016-10-31-131219 3D phase space.sagews
Views : 71
Description: Example of 3D phase space plot
# based on code from oddrobot, http://sagenb.org/home/pub/1532/
from sage.calculus.desolvers import desolve_system_rk4
x,y,t=var('x y t')
class DESolution:
def __init__(self,system,time_range,initial,stepsize=0.05,initial_points=20):
self.tvar=time_range[0]
self._times=srange(time_range[1],time_range[2],stepsize)
self.vars=[v for v,_ in initial]
self.dim=len(self.vars)
self._system=system
# check to see if we need one solution or multiple solutions
if isinstance(initial[0][1],(list,tuple)):
# multiple initial values, from the first value of each variable to the last value of each variable
start = vector([a for _,(a,b) in initial])
line = vector([b for _,(a,b) in initial])-start
self._soln = [desolve_odeint(system, ics=list(start+t*line),
times=self._times, dvars=self.vars, ivar=self.tvar) for t in srange(0,1,step=1/(initial_points-1),include_endpoint=True)]
else:
self._soln = [desolve_odeint(system, ics=[v for _,v in initial],times=self._times, dvars=self.vars, ivar=self.tvar)]

def phase_plot(self,vars=None,color='blue',**kwargs):
# find which indices the specified variables are
if vars is not None:
vars_index=[self.vars.index(v) for v in vars]
elif self.dim<=3:
vars_index=range(self.dim)
else:
vars_index=range(2)
p = Graphics()
for s in self._soln:
p+=line(s[:,vars_index],color=color,**kwargs)
# add an arrow head showing which way we are going around the phase line
half=int(s.shape[0]/2)
p+=arrow(s[half,vars_index], s[half+1,vars_index],color='red')
if len(vars_index)==2:
p.axes_labels([str(self.vars[v]) for v in vars_index])
return p

def coordinates(self,colors=None,**kwargs):
if colors is None:
colors=rainbow(len(self.vars))
p=Graphics()
legend=True
for s in self._soln:
# only want legends the first time
p+=sum(line(zip(self._times,s[:,i]), color=colors[i], legend_label=str(self.vars[i]) if legend else None,**kwargs) for i in range(self.dim))
legend=False
return p

var('x,y,z')
M = matrix([[-1,2,2],[2,2,2],[-3,-6,-6]])
show(M)
show(M.eigenvectors_right())
show(M.characteristic_polynomial())
show(M.characteristic_polynomial().roots())
show(M.eigenvalues())
F=[-x+2*y+2*z,2*x+2*y+2*z,-3*x-6*y-6*z]
solution=DESolution(F,[t,0,2],[[x,(2,-1)],(y,[2,4]),[z,[2,10]]])
solution.coordinates(['red','blue','green'])
solution.phase_plot()


(x, y, z)
$\displaystyle \left(\begin{array}{rrr} -1 & 2 & 2 \\ 2 & 2 & 2 \\ -3 & -6 & -6 \end{array}\right)$
[($\displaystyle 0$, [$\displaystyle \left(0,\,1,\,-1\right)$], $\displaystyle 1$), ($\displaystyle -2$, [$\displaystyle \left(1,\,-\frac{1}{2},\,0\right)$], $\displaystyle 1$), ($\displaystyle -3$, [$\displaystyle \left(1,\,0,\,-1\right)$], $\displaystyle 1$)]
$\displaystyle x^{3} + 5x^{2} + 6x$
[($\displaystyle 0$, $\displaystyle 1$), ($\displaystyle -2$, $\displaystyle 1$), ($\displaystyle -3$, $\displaystyle 1$)]
[$\displaystyle 0$, $\displaystyle -2$, $\displaystyle -3$]
3D rendering not yet implemented
x,y,z=var('x,y,z')

# Next we define the parameters
sigma=10
rho=28
beta=8/3

# The Lorenz equations
lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z]

# Time and initial conditions
N=250000
tmax=250
h=tmax/N
t=srange(0,tmax+h,h)
ics=[0,1,1]
sol=desolve_odeint(lorenz,ics,t,[x,y,z],rtol=1e-13,atol=1e-14)
X=sol[:,0]
Y=sol[:,1]
Z=sol[:,2]
# Plot the result
from mpl_toolkits.mplot3d import axes3d
from matplotlib import pyplot as plt
# Call the plot function if you want to plot the data
def plot():
fig = plt.figure(1)
ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)
ax.set_xlabel('X(t)')
ax.set_ylabel('Y(t)')
ax.set_zlabel('Z(t)')
plt.show()

plot()

var('a,x,y,z,w')
initial=[[x,[y,1]],[y,[z,2]],[z,[w,3]]]
p=[[x,y,z],[1,2,3]]
[v for v,_ in initial]
[a for _,(a,b) in initial]
[b for _,[_,b] in initial]
[c for _,c,_ in p]
initial[0][1]

(a, x, y, z, w) [x, y, z] [y, z, w] [1, 2, 3] [y, 2] [y, 1]