429.sagews

#Original example: Credit goes to Josh Kantor
#http://sagemath.blogspot.com/2008/01/josh-kantors-lorenz-attractor-example.html
Integer = int
RealNumber = float

def lorenz(t,y,params):
    return [params[0]*(y[1]-y[0]),y[0]*(params[1]-y[2])- y[1],y[0]*y[1]-params[2]*y[2]]

def lorenz_jac(t,y,params):
    return [ [-params[0],params[0],0],[(params[1]-y[2]),-1,-y[0]],[y[1],y[0],-params[2]],[0,0,0]]

time_range=156
T=ode_solver()
T.algorithm="bsimp"
T.function=lorenz
T.jacobian=lorenz_jac
#T.ode_solve(y_0=[.5,.5,.5],t_span=[0, time_range-1],params=[10,6.5,3],num_points=10000)
#l=[T.solution[i][1] for i in range(len(T.solution))]

#line3d(l,thickness=0.3, viewer='jmol', figsize=8)

#9.3.3
#First a plot in x,y,z
#params=[10,22,8/3]
T.ode_solve(y_0=[.5,.5,.5],t_span=[0,155],params=[30,66,8],num_points=10000)
l=[T.solution[i][1] for i in range(len(T.solution))]
#for interactive image try viewer='jmol'
line3d(l,thickness=0.3, viewer='tachyon', figsize=8)

#calculate dimensions
xs=[T.solution[i][1][0] for i in range(len(T.solution))]
ys=[T.solution[i][1][1] for i in range(len(T.solution))]
zs=[T.solution[i][1][2] for i in range(len(T.solution))]
mins=(min(xs),min(ys),min(zs))
maxs=(max(xs),max(ys),max(zs))
print([mins,maxs])

[(-48.402387945117837, -63.336232345397882, 0.40164579924863386), (50.941610214975526, 67.036089653494173, 114.31641467801353)]

xz = [(T.solution[i][1][0], T.solution[i][1][2]) for i in range(len(T.solution))]
line(xz)

xt = [((time_range-1)*i/len(T.solution), T.solution[i][1][0]) for i in range(len(T.solution))]
line(xt)

yt = [((time_range-1)*i/len(T.solution), T.solution[i][1][1]) for i in range(len(T.solution))]
line(yt)

#9.3.5
T.ode_solve(y_0=[.5,.5,.5],t_span=[0,155],params=[10,100,8/3],num_points=10000)
l=[T.solution[i][1] for i in range(len(T.solution))]
#for interactive image try viewer='jmol'
line3d(l,thickness=0.3, viewer='tachyon', figsize=8)

xz = [(T.solution[i][1][0], T.solution[i][1][2]) for i in range(len(T.solution))]
line(xz)

xt = [((time_range-1)*i/len(T.solution), T.solution[i][1][0]) for i in range(len(T.solution))]
line(xt)

yt = [((time_range-1)*i/len(T.solution), T.solution[i][1][1]) for i in range(len(T.solution))]
line(yt)