# ======================================
# Elliot Pace
# elliotlp@gmail.com
# ======================================
:
restart;
with(plots):
with(linalg):
dxdt := sigma*(y(t)-x(t)):
dydt := x(t)*(rho-z(t))-y(t):
dzdt := x(t)*y(t)-B*z(t): #using B instead of Beta because assigment of Beta through ":=" is not allowed
Lorenz:=(diff(x(t),t) = dxdt, diff(y(t),t) = dydt, diff(z(t),t) = dzdt);
fixedpoints := solve({sigma*(y-x)=0, x*(rho-z)-y=0, x*y-B*z=0}, {x,y,z}) assuming sigma>0, rho>0, B>0;
fixedpoint1 := fixedpoints[1];
fixedpoint2 := fixedpoints[2];
sigma := 10:
rho := 50:
B := 8/3:
fixedpoints := solve({sigma*(y-x)=0, x*(rho-z)-y=0, x*y-B*z=0}, {x,y,z}) assuming sigma>0, rho>0, B>0;
fixedpoint1 := fixedpoints[1];
fixedpoint2 := fixedpoints[2];
fixedpoint2 := evalf(fixedpoint2);
vfixedpoint1 := [rhs(fixedpoint1[1]), rhs(fixedpoint1[2]), rhs(fixedpoint1[3])];
vfixedpoint2 := [rhs(fixedpoint2[1]), rhs(fixedpoint2[2]), rhs(fixedpoint2[3])];
# The RootOf contains has a negative counterpart
vfixedpoint3 := [-rhs(fixedpoint2[1]), -rhs(fixedpoint2[2]), rhs(fixedpoint2[3])];
init1 := x(0)=25.0, y(0)=5.0, z(0)=25.0:
dvars := [x(t),y(t),z(t)]:
Lorenz1 := [Lorenz,init1]:
DS := dsolve(Lorenz1, numeric):
t0 := 0:
t1 := 30:
plot1 := odeplot(
DS,
dvars,
t=t0..t1,
axes=frame,
numpoints=20000,
color=blue,
thickness=1,
orientation=[20,75],
title=`Lorenz Attractor`
):
fixedpoint1plot := pointplot3d(vfixedpoint1, axes=normal, symbol=sphere, thickness=10, color=red):
fixedpoint2plot := pointplot3d(vfixedpoint2, axes=normal, symbol=sphere, thickness=10, color=green):
fixedpoint3plot := pointplot3d(vfixedpoint3, axes=normal, symbol=sphere, thickness=10, color=green):
plots[display]({plot1, fixedpoint1plot, fixedpoint2plot, fixedpoint3plot});