2018-04-13-131221 3D phase space.sagews

# 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)

$\displaystyle \left(\begin{array}{rrr} 0 & -5 & 0 \\ 5 & 0 & 0 \\ 0 & 0 & -2 \end{array}\right)$

$\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)$

$\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)$

$\displaystyle \left(\begin{array}{r} -\sqrt{2} - \frac{3}{2} \\ \sqrt{2} - \frac{3}{2} \\ -\frac{7}{2} \, \sqrt{2} \end{array}\right)$

$\displaystyle \left(\begin{array}{r} -2 \, \sqrt{2} + 2 \\ 2 \, \sqrt{2} + 2 \\ 3 \, \sqrt{2} \end{array}\right)$

$\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)$

[($\displaystyle -5 i$, [$\displaystyle \left(1,\,\frac{2}{3} i \, \sqrt{2} - \frac{1}{3},\,-\frac{1}{3} \, \sqrt{2} - \frac{2}{3} i\right)$], $\displaystyle 1$), ($\displaystyle 5 i$, [$\displaystyle \left(1,\,-\frac{2}{3} i \, \sqrt{2} - \frac{1}{3},\,-\frac{1}{3} \, \sqrt{2} + \frac{2}{3} i\right)$], $\displaystyle 1$), ($\displaystyle -2$, [$\displaystyle \left(1,\,1,\,\sqrt{2}\right)$], $\displaystyle 1$)]

$\displaystyle x^{3} + 2 x^{2} + 25 x + 50$

[($\displaystyle -5 i$, $\displaystyle 1$), ($\displaystyle 5 i$, $\displaystyle 1$), ($\displaystyle -2$, $\displaystyle 1$)]

[$\displaystyle -5 i$, $\displaystyle 5 i$, $\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]

from sage.calculus.desolvers import desolve_odeint
from sage.manifolds.utilities import set_axes_labels
x,y,z=var('x,y,z')
f=[x*(1-y),-y*(1-x),-z]
sol=desolve_odeint(f,[0.5,2,1],srange(0,10,0.1),[x,y,z])
p=Graphics()
p=line3d(zip(sol[:,0],sol[:,1],sol[:,2]))
ga = set_axes_labels(p, 'X', 'Y', 'Z', color='red')
#print(ga)
show(ga)  # the 3D frame has now axes labels
p.show()

3D rendering not yet implemented

3D rendering not yet implemented