Sharedcloud-examples / sage / 2018-04-13-131221 3D phase space.sagewsOpen in CoCalc
Author: Gerardo Emilio García Almeida
Views : 11
Description: 3D phase plane with rotation
# 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 #ADA 14 Equipo 2. var('x,y,z') A = matrix([[0,-5,0],[5,0,0],[0,0,-2]]) show(A) R = matrix([[1/2,-1/sqrt(2),1/2],[1/2,1/sqrt(2),1/2],[-1/sqrt(2),0,1/sqrt(2)]]) #R es la matriz de rotación que es composición de una rotación de -45° en el plano xz seguida de una rotación de 45° en el plano xy. show(R) show(R^-1) M=R*A*(R^-1) a=matrix([[2],[2],[-5]]) show(R*a) b=matrix([[-1],[4],[5]]) show(R*b) show(M) show(M.eigenvectors_right()) show(M.characteristic_polynomial()) show(M.characteristic_polynomial().roots()) show(M.eigenvalues()) F=[-(1/2)*x-(1/2)*(5*sqrt(2)+1)*y+(1/2)*(5-sqrt(2))*z,(1/2)*(5*sqrt(2)-1)*x-(1/2)*y-(1/2)*(sqrt(2)+5)*z,-(1/2)*(sqrt(2)+5)*x+(1/2)*(5-sqrt(2))*y-z] solution=DESolution(F,[t,0,15],[[x,(-sqrt(2)-3/2,-2*sqrt(2)+2)],(y,[sqrt(2)-3/2,2*sqrt(2)+2]),[z,[-7/sqrt(2),3*sqrt(2)]]]) solution.coordinates(['red','blue','green']) solution.phase_plot()
(x, y, z)
(050500002)\displaystyle \left(\begin{array}{rrr} 0 & -5 & 0 \\ 5 & 0 & 0 \\ 0 & 0 & -2 \end{array}\right)
(121221212122121220122)\displaystyle \left(\begin{array}{rrr} \frac{1}{2} & -\frac{1}{2} \, \sqrt{2} & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} \, \sqrt{2} & \frac{1}{2} \\ -\frac{1}{2} \, \sqrt{2} & 0 & \frac{1}{2} \, \sqrt{2} \end{array}\right)
(121212212212201212122)\displaystyle \left(\begin{array}{rrr} \frac{1}{2} & \frac{1}{2} & -\frac{1}{2} \, \sqrt{2} \\ -\frac{1}{2} \, \sqrt{2} & \frac{1}{2} \, \sqrt{2} & 0 \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \, \sqrt{2} \end{array}\right)
(232232722)\displaystyle \left(\begin{array}{r} -\sqrt{2} - \frac{3}{2} \\ \sqrt{2} - \frac{3}{2} \\ -\frac{7}{2} \, \sqrt{2} \end{array}\right)
(22+222+232)\displaystyle \left(\begin{array}{r} -2 \, \sqrt{2} + 2 \\ 2 \, \sqrt{2} + 2 \\ 3 \, \sqrt{2} \end{array}\right)
(1252212122+5252212121225212252122+521)\displaystyle \left(\begin{array}{rrr} -\frac{1}{2} & -\frac{5}{2} \, \sqrt{2} - \frac{1}{2} & -\frac{1}{2} \, \sqrt{2} + \frac{5}{2} \\ \frac{5}{2} \, \sqrt{2} - \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} \, \sqrt{2} - \frac{5}{2} \\ -\frac{1}{2} \, \sqrt{2} - \frac{5}{2} & -\frac{1}{2} \, \sqrt{2} + \frac{5}{2} & -1 \end{array}\right)
[(5i\displaystyle -5 i, [(1,23i213,13223i)\displaystyle \left(1,\,\frac{2}{3} i \, \sqrt{2} - \frac{1}{3},\,-\frac{1}{3} \, \sqrt{2} - \frac{2}{3} i\right)], 1\displaystyle 1), (5i\displaystyle 5 i, [(1,23i213,132+23i)\displaystyle \left(1,\,-\frac{2}{3} i \, \sqrt{2} - \frac{1}{3},\,-\frac{1}{3} \, \sqrt{2} + \frac{2}{3} i\right)], 1\displaystyle 1), (2\displaystyle -2, [(1,1,2)\displaystyle \left(1,\,1,\,\sqrt{2}\right)], 1\displaystyle 1)]
x3+2x2+25x+50\displaystyle x^{3} + 2 x^{2} + 25 x + 50
[(5i\displaystyle -5 i, 1\displaystyle 1), (5i\displaystyle 5 i, 1\displaystyle 1), (2\displaystyle -2, 1\displaystyle 1)]
[5i\displaystyle -5 i, 5i\displaystyle 5 i, 2\displaystyle -2]
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 = fig.add_subplot(111, projection='3d') # ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10) ax.plot(X, Y, Z) 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]