2042.sagews

x,y,z,t=var('x y z t')
a=10
b=28
c=8/3
init_x=2
init_y=3
init_z=4
p=desolve_system_rk4([a*(y-x),x*(b-z)-y,x*y-c*z],[x,y,z],ics=[0,init_x,init_y,init_z],ivar=t,end_points=20)
q1=[[i,j] for i,j,k,1 in p]
lp1=list_plot(q1,plotjoined=true)
q2=[[i,k] for i,j,k,1 in p]
lp2=list_plot(q2,plotjoined=true, color='red')
q3=[[i,1] for i,j,k,1 in p]
lp3=list_plot(q3,plotjoined=true, color='green')
lp1+lp2+lp3

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]]

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

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