Book information

Stephen Lynch
        Dynamical Systems with Applications using Mathematica®
        
This introduction to dynamical systems theory treats both continuous dynamical systems and discrete systems. Driven by numerous examples from a broad range of     disciplines and requiring only knowledge of ordinary differential equations, the text emphasizes applications and simulation utilizing the Mathematica® package.
Beginning with a tutorial guide to Mathematica®, the text thereafter is divided into three main areas. In Part I, continuous dynamical systems are considered, with examples presented from mechanical systems, chemical kinetics, electric circuits, economics, population dynamics, epidemiology, and neural networks. In Part II, both real and complex discrete dynamical systems are considered, with examples presented from population dynamics, nonlinear optics, and materials science. Part III draws on examples from both continuous and discrete systems when considering chaos control and synchronization and neural networks. Common themes such as bifurcation, bistability, chaos, fractals, instability, multistability, periodicity, and quasiperiodicity run through several chapters. Some material deals with cutting-edge published research articles and provides a useful resource for open problems in nonlinear dynamical systems.

Approximately 400 illustrations, over 400 examples, and exercises with solutions play a key role in the presentation. There are numerous simple Mathematica® programs listed throughout the text; the working notebooks may also be downloaded from the Internet at: 

						http://library.wolfram.com/infocenter/Books/Mathematics/.

Additional applications and further links of interest are also available at the author's website. 

The hands-on approach of Dynamical Systems with Applications using Mathematica® engages a wide audience of senior undergraduate and graduate students, applied mathematicians, engineers, and working scientists in various areas of the natural sciences. 

Reviews of the author’s published book Dynamical Systems with Applications using Maple®:
“The text treats a remarkable spectrum of topics…and has a little for everyone. It can serve as an introduction to many of the topics of dynamical systems, and will help even the most jaded reader, such as this reviewer, enjoy some of the interactive aspects of studying dynamics using Maple®.” - U.K. Nonlinear News
“…will provide a solid basis for both research and education in nonlinear dynamical systems.” - The Maple Reporter

ISBN 0-8176-4482-2

DYNAMICAL SYSTEMS WITH APPLICATIONS USING MATHEMATICA®
Stephen Lynch
Copyright Birkhäuser (2007)

Chapter 0: A Tutorial Introduction to Mathematica

Tutorial One: The Basics

Hold down the SHIFT key and press ENTER at the end of a line to see the output.

(*This is a comment*)

2+3-5

3 6/2-2^3

2/3+4/5

2/3+4/5 //N

Sqrt[16]

Sin[Pi]

z1=1+2I;z2=3-4I;z3=z1-z2/z1

ComplexExpand[Exp[z1]]

Factor[x^3-y^3]

Expand[%]

f=mu x(1-x) /. {mu->4,x->0.2}

Clear[mu,x]

Simplify[(x^3-y^3)/(x-y)]

Dt[x^2-3x+6,x]

D[x^3y^5,{x,2},{y,3}]

Integrate[Sin[x]*Cos[x],x]

Integrate[Exp[-x^2],{x,0,Infinity}]

Sum[1/n^2,{n,1,Infinity}]

Solve[x^2-5x+6==0,x]

Solve[x^2-5x+8==0,x]

Abs[x] /.%

Solve[{x^2+y^2==1,x+3y==0},{x,y}]

Series[Exp[x],{x,0,5}]

Limit[x/Sin[x],x->0]

f=Function[x,4 x(1-x)]

f[0.2]

fofof=Nest[f,x,3]

fnest:=NestList[f,x,4]
fofofof=fnest[[4]]

u=Table[2 i-1,{i,5}]

%^2

a={2,3,4};b={5,6,7};

3a

a.b

Cross[a,b]

Norm[a]

A={{1,2},{3,4}};B={{5,6},{7,8}};

A.B-B.A

MatrixPower[A,3]

Inverse[A]

Det[B]

Tr[B]

M={{1,2,3},{1,3,2},{3,2,1}}

Eigenvalues[M]

Eigenvectors[M]

LaplaceTransform[t^3,t,s]

InverseLaplaceTransform[6/s^4,s,t]

FourierTransform[t^4 Exp[-t^2],t,w]

InverseFourierTransform[%,w,t]

Quit[ ];

Tutorial Two: Plots and Differential Equations

(* Tutorial 2*)

Plot[Sin[x],{x,-Pi,Pi}]

Plot[{Cos[x],Exp[-.1x] Cos[x]},{x,0,60}]

Plot[Exp[-.1t] Sin[t],{t,0,60},AxesLabel->{"t","Current"},PlotRange->{-1,1}]

ContourPlot[y^2/2-x^2/2+x^4/4,{x,-2,2},{y,-2,2}]

ContourPlot[y^2/2-x^2/2+x^4/4,{x,-2,2},{y,-2,2},ContourShading->False]

ContourPlot[y^2/2-x^2/2+x^4/4,{x,-2,2},{y,-2,2},Contours->20]

Plot3D[y^2/2-x^2/2+x^4/4,{x,-2,2},{y,-2,2}]

ParametricPlot[{t^3-4t,t^2},{t,-3,3}]

ParametricPlot3D[{Sin[t],Cos[t],t/3},{t,-10,10}]

<<Graphics`ImplicitPlot`

ImplicitPlot[2x^2+3y^2==12,{x,-3,3},{y,-3,3}]

DSolve[x'[t]==-x[t]/t,x[t],t]

DSolve[{x'[t]==-t/x[t],x[0]==1},x[t],t]

DSolve[{x''[t]+5x'[t]+6x[t]==10Sin[t],x[0]==0,x'[0]==0},x[t],t]

DSolve[{x'[t]==3x[t]+4y[t],y'[t]==-4x[t]+3y[t]},{x[t],y[t]},t]

DSolve[{x'[t]==x[t],y'[t]==y[t],z'[t]==-z[t]},{x[t],y[t],z[t]},t]

NDSolve[{x'[t]==x[t](.1-.01x[t]),x[0]==50},x,{t,0,100}]

Plot[Evaluate[x[t]/.%],{t,0,100},PlotRange->All]

v=NDSolve[{x'[t]==.1x[t]+y[t],y'[t]==-x[t]+.1y[t],x[0]==.01,y[0]==0},{x[t],y[t]},{t,0,50}]

ParametricPlot[{x[t],y[t]}/.v,{t,0,50} , PlotRange->All, PlotPoints->1000]

w=NDSolve[{x'[t]==z[t]-x[t],y'[t]==-y[t],z'[t]==z[t]-17x[t]+16,x[0]==.8,y[0]==.8,z[0]==.8},{x[t],y[t],z[t]},{t,0,20}]

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.w],{t,0,20} , PlotPoints->1000,PlotRange->All]

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
vanderpol=NDSolve[{Derivative[1][x][t]==y[t],Derivative[1][y][t]==1000*(1-x[t]^2)*y[t]-x[t],x[0]==2,y[0]==0},{x,y},{t,5000}];
T=First[InterpolatingFunctionCoordinates[First[x/.vanderpol]]];
ListPlot[Transpose[{x[T],y[T]}/.First[vanderpol]],PlotRange->All]

Quit[ ];

Tutorial Three: The Manipulate Command and Simple Mathematica Programs

The Manipulate Command: To see animation commands click on the plus icon at the end of the slider.

Example 1: Solving simultaneous equations as the parameter a varies. (See Exercise 6 in Chapter 3).

Manipulate[Solve[{x (1-y-a x)==0,y (-1+x-a y)==0}],{a,0,1}]

Showing the solution lying wholly in the first quadrant.

Manipulate[Plot[{1-a x,(x-1)/a},{x,0,2}],{a,0.01,1}]

Example 2: Evaluating functions of functions to the fifth iteration.

f=Function[x,4 x (1-x)];

Manipulate[Expand[Nest[f,x,d]],{d,1,5,1}]

Example 3: Solving a differential equation as two parameters vary. In this case, the damping and restoring terms when modeling a simple pendulum.

Manipulate[DSolve[{x''[t]+b x'[t]+c x[t]==0,x[0]==1,x'[0]==0},x[t],t],{b,0,2},{c,0,5}]

The solution curves.

Manipulate[soln=NDSolve[{x''[t]+b x'[t]+c x[t]==0,x[0]==1,x'[0]==0},x[t],{t,0,100}];
Plot[Evaluate[x[t]/.soln],{t,0,100},PlotRange->All],{b,0,5},{c,0,5}]

Simple Mathematica Programs

Modules and local variables

Norm3d[a_,b_,c_]=Module[{t},t=Sqrt[a^2+b^2+c^2]]

Norm3d[3,4,5]

The Do command

F[1]=1;F[2]=1;Nmax=10;
Do[F[i]=F[i-1]+F[i-2],{i,3,Nmax}]

Table[F[i],{i,1,Nmax}]

If, then, else construct

Num=1000;If[Num>0,Print["Num is positive"],If[Num==0,Print["Num is zero"],Print["Num is negative"]]]

Conditional statements with more than two alternatives

r[x_]=Switch[Mod[x,3],0,a,1,b,2,c]

r[7]

Which command

T[x_]=Which[0≤x<1/2,mu x,1/2≤x≤1,mu(1-x)]       (*Tent function*)

T[4/5]

For loops

For[i=1;t=x,i<4,i++,t=mu t (1-t);Print[Factor[t]]]         (*f^n[x]*)

While command

a=12348;b=14238;

While[b≠0,d=Mod[a,b];a=b;b=d]

a

Quit[ ];

Plotting solution curves

ode1[x0_,y0_]:=NDSolve[{x'[t]==y[t] ,y'[t]==-25 x[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,0,6}];
sol1[1]=ode1[1,0] ;
ode2[x0_,y0_]:=NDSolve[{x'[t]==y[t] ,y'[t]==-2 y[t]-25 x[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,0,6}];
sol2[1]=ode2[1,0] ;
p1=Plot[Evaluate[Table[{x[t]}/.sol1[i],{i,1}]],{t,0,6}, PlotRange->All, PlotPoints->100,PlotStyle->Dashing[{.02}]];
p2=Plot[Evaluate[Table[{x[t]}/.sol2[i],{i,1}]],{t,0,6}, PlotRange->All, PlotPoints->100];
Show[{p1,p2},PlotRange->All,AxesLabel->{"t","x"},Axes->True,TextStyle->{FontSize->15}]

Quit[ ];

Chapter 1: Differential Equations

Example 1

Solving a simple separable differential equation(dx/dt)=-(t/x).

DSolve[x'[t]==-t/x[t],x[t],t]

Example 8 (Figure 1.6)

Plotting the solution curve to an initial value problem (dc/dt)=k((((Subscript[a, 0])-c))^2)((((Subscript[b, 0])-(c/2)))^2), c(0)=0.

k=.00713;a0=4;b0=1;sol=NDSolve[{c'[t]==k (a0-c[t]^2)(b0-c[t]/2),c[0]==0},c,{t,0,700}]

Plot[Evaluate[c[t]/.sol],{t,0,700}]

Example 10

Solving an initial value problem (((d^2)I)/(dt^2))+5(dI/dt)+6I=10sin(t), I(0)=I'(0)=0.

DSolve[{i''[t]+5i'[t]+6i[t]==10Sin[t],i[0]==0,i'[0]==0},i[t],t]

Exercise 7

Solving a system of three differential equations (dx/dt)=-ax, (dy/dt) =ax-by, (dz/dt)=by, x(0)=M, y(0)=z(0)=0.

DSolve[{x'[t]==-a x[t],y'[t]==a x[t]-b y[t],z'[t]==b y[t],x[0]==M,y[0]==0,z[0]==0},{x[t],y[t],z[t]},t];

Simplify[%]

Exercise 10(b)

Solving an initial value problem (dx/dt)=2(Sqrt[x]), x((Subscript[t, 0]))=(Subscript[x, 0]).

DSolve[{x'[t]==2 Sqrt[x[t]],x[t0]==x0},x[t],t]
Simplify[%]

Quit[ ];

Chapter 2: Planar Systems

Example 4 (Figure 2.8)

Plotting simple phase plane portraits (dx/dt)=2x+y, (dy/dt) =x+2y. Eight trajectories are plotted with a vector field superimposed.

<<Graphics`PlotField`

p1=PlotVectorField[{2x+y,x+2y},{x,-3,3},{y,-3,3}] ;

ode1[x0_,y0_]:=NDSolve[{x'[t]==2x[t]+y[t],y'[t]==x[t]+2y[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,-3,3}];

sol[1]=ode1[1,1];sol[2]=ode1[1,-1];sol[3]=ode1[-1,-1];sol[4]=ode1[-1,1];sol[5]=ode1[3,1];sol[6]=ode1[1,3];sol[7]=ode1[-1,-3];sol[8]=ode1[-3,-1];

p2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,8}]],{t,-3,3}, PlotRange->{{-3,3},{-3,3}},PlotPoints->100,AxesLabel->{"x","y"}];
Show[{p1,p2},PlotRange->{{-3,3},{-3,3}},AxesLabel->{"x","y"},Axes->True]

Example 8

Locating the critical points of the system (dx/dt)=y, (dy/dt) =x(1-(x^2))+y.

Solve[{y==0,x (1-x^2)+y==0},{x,y}]

Example 8 (Figure 2.12)

Plotting the phase portrait of (dx/dt)=y, (dy/dt) =x(1-(x^2))+y.

<<Graphics`PlotField`

p1=PlotVectorField[{y,x (1-x^2)+y},{x,-3,3},{y,-3,3}] ;

ode1[x0_,y0_]:=NDSolve[{x'[t]==y[t],y'[t]==x[t] (1-(x[t])^2)+y[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,-10,10}];

sol[1]=ode1[0.1,0.16] ;sol[2]=ode1[0.1,-0.06];sol[3]=ode1[-0.1,-0.16];
sol[4]=ode1[-0.1,0.06];sol[5]=ode1[1.3,0] ;                     sol[6]=ode1[-1.3,0];

p2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,6}]],{t,-10,10}, PlotRange->{{-3,3},{-3,3}},PlotPoints->100,AxesLabel->{"x","y"}];

Show[{p1,p2},PlotRange->{{-3,3},{-3,3}},AxesLabel->{"x","y"},Axes->True]

Quit[ ];

Chapter 3: Interacting Species

Example 3 (Figure 3.5(a))

The Holling-Tanner model: (dx/dt)=x(1-(x/7))-((6xy)/(7+7x)), (dy/dt) =0.2y(1-((0.5y)/x)).

<<Graphics`PlotField`

p1=PlotVectorField[{x (1-x/7)-6x y/(7+7x),.2 y (1-.5 y/x)},{x,.1,7},{y,.1,5},DisplayFunction->Identity] ;

ode1[x0_,y0_]:=NDSolve[{x'[t]==x[t] (1-x[t]/7)-6 x[t] y[t]/(7+7 x[t]),y'[t]==.2 y[t] (1-.5 y[t]/x[t]),x[0]==x0,y[0]==y0},{x[t],y[t]},{t,0,200}];
sol[1]=ode1[7,.1] ;sol[2]=ode1[1.1,2];

p2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,2}]],{t,0,200}, PlotRange->{{0,7},{0,5}}, PlotPoints->100,AxesLabel->{"x","y"},DisplayFunction->Identity];

Show[{p1,p2},PlotRange->{{0,7},{0,5}},AxesLabel->{"x","y"},Axes->True,DisplayFunction->$DisplayFunction]

Example 3
           Time series of the Holling-Tanner model.

p3=Plot[Evaluate[Table[{x[t]}/.sol[i],{i,1}]],{t,0,200}, PlotRange->All, PlotPoints->100,PlotStyle->Dashing[{.02}],DisplayFunction->Identity];

p4=Plot[Evaluate[Table[{y[t]}/.sol[i],{i,1}]],{t,0,200}, PlotRange->All, PlotPoints->100,DisplayFunction->Identity];

Show[{p3,p4},PlotRange->All,AxesLabel->{"t","x,y"},Axes->True,DisplayFunction->$DisplayFunction]

Quit[ ];

Chapter 4: Limit Cycles

The Fitzhugh-Nagumo Oscillator (Figure 4.3)

Phase portrait and isoclines: (du/dt)=-u(u-θ)(u-1)-v+w,(dv/dt) =ϵ(u-γv).

theta=.14;gamma=2.54;epsilon=.01;omega=.112;

p1=Plot[{-u (u-theta) (u-1)+omega,u/gamma},{u,-.4,1.2},DisplayFunction->Identity];

sol=NDSolve[{u'[t]==-u[t] (u[t]-theta) (u[t]-1)-v[t]+omega,v'[t]==epsilon (u[t]-gamma v[t]),u[0]==.5,v[0]==.1},{u[t],v[t]},{t,50,150}];

p2=ParametricPlot[Evaluate[Table[{u[t],v[t]}/.sol]],{t,50,150}, PlotRange->{{-.4,1.2},{0,.4}}, PlotPoints->100,AxesLabel->{"u","v"}];
Show[{p1,p2},PlotRange->{{-.4,1.2},{0,.4}},AxesLabel->{"u","v"},Axes->True]

Example 8. The Lidstedt-Poincaré technique.

epsilon="epsilon";

SetAttributes[{w1,epsilon},Constant]

x=x0+epsilon*x1+epsilon^2*x2;

Collect[(1+2 epsilon w1+epsilon^2 (w1^2+2 w2))Dt[x,{t,2}]+x-epsilon x^3,epsilon]

DSolve[{x0''[t]+x0[t]==0,x0[0]==1,x0'[0]==0},x0[t],t]

DSolve[{x1''[t]+x1[t]==Cos[t]^3+2 w1 Cos[t],x1[0]==0,x1'[0]==0},x1[t],t]

Simplify[%]

Quit[ ];

Chapter 5: Hamiltonian Systems, Lyapunov Functions, and Stability

Example 6 (Figure 5.8(a)).

Surface plot.

gama=.7;
Plot3D[-(x^2+y^2)-(4/(gama Pi^2))(Log[Cos[Pi x/2]]+Log[Cos[Pi y/2]]),{x,-1,1},{y,-1,1},ViewPoint->{1,2,1},PlotPoints->50,AxesLabel->{"(Subscript[a, 1])","(Subscript[a, 2])","V"}]

Example 6 (Figure 5.8(b)).

Contour plot.

ContourPlot[- ((x^2+y^2)+(4/(gama Pi^2))(Log[Cos[Pi x/2]]+Log[Cos[Pi y/2]])),{x,-1,1},{y,-1,1},PlotPoints->50,ContourShading->False,Contours->50,Axes->True,AxesLabel->{"(Subscript[a, 1])","(Subscript[a, 2])"}]

Example 6 (Figure 5.8(c)).

Density plot.

DensityPlot[-(x^2+y^2)-(4/(gama Pi^2))(Log[Cos[Pi x/2]]+Log[Cos[Pi y/2]]),{x,-1,1},{y,-1,1},PlotPoints->50,ColorFunction->(GrayLevel[#^.4]&),Mesh->False,Axes->True,AxesLabel->{"(Subscript[a, 1])","(Subscript[a, 2])"}]

Quit[ ];

Chapter 6: Bifurcation Theory

Example 2

Taylor series expansion.

Series[1/(1+2 u),{u,0,3}]

Exercise 8(b)

Animation of a simple curve as a parameter varies.

Manipulate[Plot[r ((r-1)^2-mu r),{r,-3,3},PlotRange->{-3,3}],{mu,-1,1}]

Manipulate[sol=NDSolve[{x'[t]==y[t]+mu x[t]-x[t] y[t]^2,y'[t]==mu y[t]-x[t]-y[t]^3,x[0]==.5,y[0]==.5},{x[t],y[t]},{t,0,80}];ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol]],{t,0,80}, PlotRange->{{-1,1},{-1,1}}, PlotPoints->100],{mu,-1,0.3}]

Quit[ ];

Chapter 7: Three-Dimensional Autonomous Systems and Chaos

The Lorenz Attractor (Figure 7.11)

The Lorenz equations: (dx/dt)=10(y-x), (dy/dt) =-xz+28x-y, (dz/dt)=xy-(8/3)z.

Lorenz=NDSolve[{x'[t]==10 (y[t]-x[t]),y'[t]==-x[t] z[t]+28 x[t]-y[t],z'[t]==x[t] y[t]-(8/3) z[t],x[0]==15,y[0]==20,z[0]==30},{x,y,z},{t,0,100},MaxSteps->Infinity];

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.Lorenz],{t,0,100},PlotPoints->10000]

The Double Scroll Attractor (Figure 7.14)

Chua's equations: (dx/dt)=15(y-x-(cx+(1/2)(d-c)(|x+1|-|x-1|))), (dy/dt) =x-y+z , (dz/dt)=-by.

Chua=NDSolve[{x'[t]==15 (y[t]-x[t]-((-5/7) x[t]+(1/2) (-(8/7)-(-5/7)) (Abs[x[t]+1]-Abs[x[t]-1]))),y'[t]==x[t]-y[t]+z[t],z'[t]==-25.58 y[t],x[0]==-1.6,y[0]==0,z[0]==1.6},{x,y,z},{t,0,100},MaxSteps->Infinity];

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.Chua],{t,0,100},PlotPoints->10000]

Exercise 6

The Chapman Cycle: A stiff system of differential equations used to model the production of ozone.

Chapman=NDSolve[{x'[t]==2 3 10^(-12) y[t]+5.5 10^(-4) z[t]-1.22 10^(-33) x[t] y[t] 9 10^17-6.86 10^(-16) x[t] z[t],y'[t]==5.5 10^(-4) z[t]+2 6.86 10^(-16) x[t] z[t]-3 10^(-12) y[t]-1.22 10^(-33) x[t] y[t] 9 10^(17),z'[t]==1.22 10^(-33) x[t] y[t] 9 10^17-5.5 10^(-4) z[t]-6.86 10^(-16) x[t] z[t],x[0]==4.6806 10^7,y[0]==6.999 10^16,z[0]==6.5396 10^12},{x,y,z},{t,0,10^8},MaxSteps->Infinity];

Singlet=x[10^8] /. Chapman

Oxygen=y[10^8] /. Chapman

Ozone=z[10^8] /. Chapman

Quit[ ];

Chapter 8: Poincaré Maps and Nonautonomous Systems in the Plane

Example 1

Poincaré first returns for(dr/dt)=-(r^2).

returns=NDSolve[{r'[t]==-(r[t])^2,r[0]==1},r,{t,0,16 Pi}]

Table[Evaluate[r[2i Pi]],{i,8}] /. returns

Example 5 (Figure 8.5(e))

Hamiltonian system: ((Subscript[dp, 1])/dt)=-(Sqrt[2])(Subscript[q, 1]), ((Subscript[dq, 1])/dt) =(Sqrt[2])(Subscript[p, 1]), ((Subscript[dp, 2])/dt)=-(Subscript[q, 2]), ((Subscript[dq, 2])/dt)=-(Subscript[p, 2]).

Ham1=NDSolve[{p1'[t]==-Sqrt[2] q1[t],q1'[t]==Sqrt[2] p1[t],p2'[t]==-q2[t],q2'[t]==p2[t],p1[0]==.5,p2[0]==1.5,q1[0]==.5,q2[0]==0},{p1,p2,q1,q2},{t,0,1000},MaxSteps->Infinity];

ParametricPlot3D[Evaluate[{p1[t],q1[t],q2[t]} /.Ham1],{t,0,200},AxesLabel->{"p1","q1","p2"},PlotPoints->10000]

Example 5 (Figure 8.5(f))

Two-dimensional surface of section.

ParametricPlot3D[Evaluate[{p1[t],q1[t],q2[t]} /.Ham1],{t,0,1000},AxesLabel->{"p1","q1"," "},PlotRange->{{-1,1},{-1,1},{-0.01,0.01}},PlotPoints->10000]

The Duffing System (Figure 8.11(b))

A nonautonomous system: (dx/dt)=y, (dy/dt)=x-ky-(x^3)+Acos(ωt).

k=0.3;A=0.5;omega=1.25;

Duffing=NDSolve[{x'[t]==y[t],y'[t]==x[t]-k y[t]-(x[t])^3+A Cos[omega t],x[0]==1,y[0]==0},{x[t],y[t]},{t,0,200},MaxSteps->Infinity]

ParametricPlot[Evaluate[{x[t],y[t]} /.Duffing],{t,0,200},PlotPoints->10000]

The Duffing System (Figure 8.11(b))

Poincaré section.

i=4000;
Poincare=NDSolve[{x'[t]==y[t],y'[t]==x[t]-k y[t]-(x[t])^3+A Cos[omega t],x[0]==1,y[0]==0},{x[t],y[t]},{t,0,2 i  Pi/omega},MaxSteps->Infinity]

ListPlot[Table[Flatten[{x[t],y[t]} /. Poincare],{t,0,2 i Pi/omega,2 Pi/omega}],PlotStyle->PointSize[0.001],AxesLabel->{"x","y"}]

Exercise 10(a) (Figure 19.2(a))

Bifurcation diagram for a periodically forced pendulum.

k=0.1;omega=1.25;
M=1200;step=0.0001;MaxAmp=M*step;a=1;b=0;
For[i=1,i<=M,i++,A=step*i;pp=NDSolve[{x'[t]==y[t],y'[t]==x[t]-k y[t]-(x[t])^3+A Cos[omega t],x[0]==a,y[0]==b},{x,y},{t,0,2 Pi/omega},MaxSteps->100000];a=Flatten[x[2 Pi/omega] /. pp];b=Flatten[y[2 Pi/omega] /. pp];rampup[i]=Sqrt[a^2+b^2]]
p1=ListPlot[Table[Flatten[{m*step,rampup[m]}],{m,1,M}]];
For[i=1,i<=M,i++,A=MaxAmp-step*i;pp=NDSolve[{x'[t]==y[t],y'[t]==x[t]-k y[t]-(x[t])^3+A Cos[omega t],x[0]==a,y[0]==b},{x,y},{t,0,2 Pi/omega},MaxSteps->100000];a=Flatten[x[2 Pi/omega] /. pp];b=Flatten[y[2 Pi/omega] /. pp];rampdown[i]=Sqrt[a^2+b^2]]
p2=ListPlot[Table[Flatten[{MaxAmp-m*step,rampdown[m]}],{m,1,M}]];
Show[{p1,p2},PlotRange->All,AxesLabel->{"Γ","r"},AxesOrigin->{0,0.4},TextStyle->{FontSize->15},Ticks->{{0,0.04,0.08,0.12},Automatic}]

Quit[ ];

Chapter 9: Local and Global Bifurcations

Section 9.1

Computation of focal values when (dx/dt)=y-(Subscript[a, 1])x-(Subscript[a, 2])(x^2)-(Subscript[a, 3])(x^3)-(Subscript[a, 4])(x^4), (dy/dt) =-x-(Subscript[b, 2])(x^2)-(Subscript[b, 3])(x^3).

v[1,0]=0;v[0,1]=0;v[2,0]=1/2;v[1,1]=0;v[0,2]=1/2;v[0,0]=0;

vv=Sum[Sum[v[i,j-i] x^i y^(j-i),{i,0,j}],{j,0,20}];

xdot=y-a2 x^2-a3 x^3-a4 x^4;

dydt=-x-b2 x^2-b3 x^3;

dd=Expand[D[vv,x] xdot+D[vv,y] dydt-eta4 (x^2+y^2)^2-eta6(x^2+y^2)^3];

Do[eq[i]=Coefficient[Coefficient[dd,x,i],y,3-i],{i,0,3}];

SolsV3=Solve[{eq[0]==0,eq[1]==0,eq[2]==0,eq[3]==0},{v[3,0],v[2,1],v[1,2],v[0,3]}]

Do[eq[i]=Coefficient[Coefficient[dd,x,i],y,4-i],{i,0,4}];

SolsV4=Solve[{eq[0]==0,eq[1]==0,eq[2]==0,eq[3]==0,eq[4]==0,v[2,2]==0},{v[4,0],v[3,1],v[2,2],v[1,3],v[0,4],eta4}] /.SolsV3

a3=2 a2 b2/3;

Do[eq[i]=Coefficient[Coefficient[dd,x,i],y,5-i],{i,0,5}];

SolsV5=Simplify[Solve[{eq[0]==0,eq[1]==0,eq[2]==0,eq[3]==0,eq[4]==0,eq[5]==0},{v[5,0],v[4,1],v[3,2],v[2,3],v[1,4],v[0,5]}] /.SolsV4 /.SolsV3]

Do[eq[i]=Coefficient[Coefficient[dd,x,i],y,6-i],{i,0,6}];

SolsV6=Simplify[Solve[{eq[0]==0,eq[1]==0,eq[2]==0,eq[3]==0,eq[4]==0,eq[5]==0,eq[6]==0,v[4,2]+v[2,4]==0},{v[6,0],v[5,1],v[4,2],v[3,3],v[2,4],v[1,5],v[0,6],eta6}] /.SolsV5 /.SolsV4/.SolsV3];

Flatten[eta6 /.SolsV6]

Quit[];

Example 4

Multivariate polynomial division.

PolynomialReduce[x^4+y^4+z^4,{x^2+y,z^2 y-1,y-z^2},{x,y,z}]

PolynomialReduce[x^4+y^4+z^4,{y-z^2,z^2 y-1,x^2+y},{x,y,z}]

Example 6

Computing Groebner bases.

polys={y^2-x^3+x,2 x^2 y^2+x^4+4 x^3+y^4-12 x y^2};
GroebnerBasis[polys,{x,y},MonomialOrder->Lexicographic]

GroebnerBasis[polys,{x,y},MonomialOrder->DegreeLexicographic]

Example 6

Solving multivariate polynomial equations.

N[Solve[{x+y^2-x^3==0,2 x^2 y^2+x^4+4 x^3+y^4-12 x y^2==0},{x,y}],4]

Example 7

Computing a Groebner basis.

polys={-a1,2 a2 b2-3 a3,5b2(2a4-b3 a2),-5b2(92b2^2a4-99b3^2a2+1520a2^2a4-760a2^3b3-46b2^2b3 a2+198b3 a4),-b2(14546b2^4a4+105639a2^3b3^2+96664a2^3b2^2b3-193328a2^2b2^2a4-891034a2^4a4+445517a2^5b3+211632a2 a4^2-317094a2^2b3 a4-44190b2^2b3 a4+22095b2^2b3^2a2-7273b2^4b3 a2+5319b3^3a2-10638b3^2a4) };
GroebnerBasis[polys,{a1,a2,a3,a4,b2,b3},MonomialOrder->Lexicographic]

Quit[ ];

Chapter 10: The Second Part of Hilbert's Sixteenth Problem

A  Liénard System with a Large Parameter (Figure 10.12(a))

p1=Plot[-x+x^3,{x,-2,2}];

sol=NDSolve[{x'[t]==10 y[t]-10 (-x[t]+x[t]^3),y'[t]==-x[t]/10,x[0]==1,y[0]==0},{x[t],y[t]},{t,50,100}];

p2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol]],{t,50,100}, PlotRange->{{-2,2},{-1,1}}, PlotPoints->100,AxesLabel->{"x","y"}];
Show[{p1,p2},PlotRange->{{-2,2},{-1,1}},AxesLabel->{"x","y"},Axes->True]

Quit[ ];

Chapter 11: Linear Discrete Dynamical Systems

Example 1

Solving recurrence relations.

Savings=RSolve[{x[n+1]==(1+(3/100)) x[n],x[0]==10000},x[n],n]

n=5;

N[Savings,7]

Clear[n]

Example 2(i)

Solving a second-order linear difference equation.

RSolve[{x[n+2]==x[n+1]+6 x[n],x[0]==1,x[1]==2},x[n],n]

Example 4

Leslie  Matrices

L={{0,3,1},{.3,0,0},{0,.5,0}}

X0={{1000},{2000},{3000}}

X10=MatrixPower[L,10].X0

Eigenvalues[L]

Eigenvectors[L]

Quit[ ];

Chapter 12: Nonlinear Discrete Dynamical Systems

Example 1

Iteration of the tent map.

mu=2;
T[x_]=Which[0≤x≤1/2,mu x,1/2≤x≤1,mu (1-x)];

NestList[T,1/11,20]
(* See Table 2.1.*)

NestList[T,0.2,20]

NestList[T,.2001,20]

Iterations=20;x0=1/11;

NL=NestList[T,x0,Iterations]

Example 3 (Figure 12.6(d))

Graphical iteration of the tent map.

x=Drop[Prepend[Flatten[Partition[NL,2,1]],x0],-Iterations];
y=Drop[Prepend[Drop[Flatten[Partition[NL,2,1]],1],0],-(Iterations-1)];

p1=ListPlot[Table[{x[[i]],y[[i]]},{i,1,Iterations+1}],PlotJoined->True,PlotStyle->RGBColor[1,0,0]];

p2=Plot[{x,T[x]},{x,0,1}];

Show[{p1,p2},PlotRange->All,AxesLabel->{"x","T"},Axes->True,AspectRatio->Automatic]

Computing the Lyapunov Exponent of the Logistic Map

Clear[mu,x,i];mu=4;

fmu[x_]=N[mu x (1-x)];

Itermax=50000;

p=NestList[fmu,.1,Itermax];

pp=Log[Abs[mu (1-2 p)]];

Total[pp/(Itermax+1)]

Plotting a Bifurction Diagram of the Logistic Map using the First Iterative Method Figure (12.15)

x0=0.4;stepsize=0.001;

Itermax=200;

iterate=Compile[{mu},Map[{mu,#}&,Union[Drop[NestList[mu # (1-#)&,x0,Itermax],100]]]];

pts=Flatten[Table[iterate[mu],{mu,0,4,stepsize}],1];

ListPlot[pts,PlotStyle->{PointSize[0.001],RGBColor[1,0,0]},AxesLabel->{"μ","x"}]

Quit[ ];

Iteration of the Hénon Map (12.22(b))

a=1.2;b=0.4;Itermax=10000;

Henon[{x_,y_}]={1+y-a x^2,b x};

HenonMap=Drop[NestList[Henon,{0.1,0},Itermax],100];

ListPlot[HenonMap,AxesLabel->{"x","y"},Axes->None,PlotStyle->{PointSize[0.001],RGBColor[0,1,0]},Frame->True]

Quit[ ];

Chapter 13: Complex Iterative Maps

A Julia Set (Figure 13.1(c))

Coordinates={Complex[x_,y_]->{x,y}};

JuliaSet[c_,z0_,Npoints_]:=Module[{Nc=N[c],Nz0=N[z0]},NestList[(1-2Random[Integer]) Sqrt[#-Nc]&,Nz0,Npoints]/.Coordinates];
JuliaData=Drop[Evaluate[JuliaSet[-1+0 I,0+5I,2^18]],100];
ListPlot[JuliaData,Axes->False,AspectRatio->Automatic,PlotStyle->PointSize[0.001]]

Quit[ ];

The Color Mandelbrot Set (Figure 13.2)

xmin=-2.2;xmax=0.6;ymin=-1.2;ymax=1.2;Itermax=50;stepsize=0.005;
Counter =Compile[{x,y,{Itermax, _Integer}}, Module[{w,z,count=0},w=x + I y;z = w;While[(Abs[z]<4.0) && (count≤Itermax),++count;z = z^2+ w;];Return[count]] ];
MandelbrotData = Table[Counter[x,y,Itermax],{y,ymin,ymax,stepsize},{x,xmin,xmax,stepsize}];

ListDensityPlot[MandelbrotData,Frame->False,Mesh->False,MeshRange->{{xmin,xmax},{ymin,ymax}},AspectRatio->Automatic,ColorFunction->(Hue[Log[1+#]]&)]

Quit[ ];

Chapter 14: Electromagnetic Waves and Optical Resonators

Chaotic Attractor of the Ikeda Map (Figure 14.11(b))

A=10;B=0.15;Itermax=10000;
Ikeda[{x_,y_}]={A+B (x Cos[x^2+y^2]-y Sin[x^2+y^2]),B (x Sin[x^2+y^2]+y Cos[x^2+y^2])};

IkedaMap=Drop[NestList[Ikeda,{A,0},Itermax],100];

ListPlot[IkedaMap,AxesLabel->{"Re(E_n)","Im(E_n)"},Axes->None,PlotStyle->{PointSize[0.001],RGBColor[0,0,1]},Frame->True,AspectRatio->Automatic]

Bifurcation Diagram for the Simple Fiber Ring Resonator using the Second Iterative Method

MaxPower=120;kappa=0.0225;Const=0.345913;Initial=0+0 I;
HalfN=19999;phi=0;N1=1+HalfN;E1[1]=0;Itermax=2 HalfN+1;Esqr[1]=0;
For[n=1,n≤HalfN,n++,{E2[n+1]=E1[n] Exp[I (Abs[Const E1[n]]^2-phi)],E1[n+1]=I Sqrt[1-kappa] Sqrt[n MaxPower/N1]+Sqrt[kappa] E2[n+1],Esqr[n+1]=Abs[E1[n+1]]^2}]
p1=ListPlot[Table[Flatten[{n*MaxPower/N1,Esqr[n]}],{n,1,HalfN}],PlotStyle->PointSize[0.001]];
For[n=N1,n≤Itermax,n++,{E2[n+1]=E1[n] Exp[I (Abs[Const E1[n]]^2-phi)],E1[n+1]=I Sqrt[1-kappa] Sqrt[2 MaxPower-n MaxPower/N1]+Sqrt[kappa] E2[n+1],Esqr[n+1]=Abs[E1[n+1]]^2}];
p2=ListPlot[Table[Flatten[{n*MaxPower/N1,Esqr[Itermax+1-n]}],{n,1,HalfN}],PlotStyle->PointSize[0.001]];
Show[{p1,p2},PlotRange->All,AxesLabel->{"Input","Output"},TextStyle->{FontSize->15},AspectRatio->Automatic]

Quit[ ];

Chapter 15: Fractals and Multifractals

The Koch Curve (Figure 15.2)

k=6;Itermax=4^k;h=3^(-k);

x[1]=0;y[1]=0;

angle[1]=0;angle[2]=Pi/3;angle[3]=-Pi/3;angle[4]=0;

For[a=1,a≤Itermax,a++,m=a-1;ang=0;For[b=1,b≤k,b++,segment[b]=Mod[m,4];m=Floor[m/4];r=segment[b]+1;ang=ang+angle[r]];x[a+1]=x[a]+h Cos[ang];y[a+1]=y[a]+h Sin[ang]]

ListPlot[Table[{x[i],y[i]},{i,1,Itermax+1}],PlotJoined->True,Axes->False,AspectRatio->Automatic]

The Sierpiński Triangle (Figure 15.6)

V1={0,0};V2={4,0};V3={2,2 Sqrt[3]};

Itermax=100000;Half=1/2;P[1]={0,0};

P[n_]=For[n=1,n≤Itermax,n++,r=Random[ ];If[0≤r≤1/3,P[n+1]=N[P[n]+(V1-P[n]) Half],If[1/3<r≤2/3,P[n+1]=N[P[n]+(V2-P[n]) Half],P[n+1]=N[P[n]+(V3-P[n]) Half]]]]

Sierpinski=Table[P[i],{i,Itermax}];

ListPlot[Sierpinski,Axes->False,AspectRatio->Automatic,PlotStyle->{PointSize[0.001],RGBColor[1,0,0]}]

Quit[];

Barnsley's Fern (Figure 15.7)

Itermax=100000;P[1]={0.5,0.5};

T[{x_,y_,a_,b_,c_,d_,e_,f_}]={a*x+b*y+c,d*x+e*y+f};

P[n_]=For[n=1,n≤Itermax,n++,r=Random[ ];If[0≤r<0.05,P[n+1]=T[{P[n][[1]],P[n][[2]],0,0,0,0,0.2,0}],If[r<0.86,P[n+1]=T[{P[n][[1]],P[n][[2]],0.85,0.05,0,-0.04,0.85,1.6}],
If[r<0.93,P[n+1]=T[{P[n][[1]],P[n][[2]],0.2,-0.26,0,0.23,0.22,1.6}],
P[n+1]=T[{P[n][[1]],P[n][[2]],-0.15,0.28,0,0.26,0.24,0.44}]]]]]

BarnsleyFern=Table[P[i],{i,Itermax}];

ListPlot[BarnsleyFern,Axes->False,AspectRatio->Automatic,PlotStyle->{PointSize[0.001],RGBColor[0,1,0]}]

Multifractal τ Curve (Figure 15.16(a))

Plot[Log[(1/9)^x+(8/9)^x]/Log[3],{x,-20,20},AxesLabel->{"q","tau(q)"}]

Multifractal (Subscript[D, q]) Curve (Figure 15.16(b))

Plot[(Log[(1/9)^x+(8/9)^x]/Log[3]) 1/(1-x),{x,-20,20}]

Multifractal f(α) Curve (Figure 15.16(c))

k=500;

ParametricPlot[{(s Log[1/9]+(k-s) Log[8/9])/(Log[3^(-k)]),-Log[Binomial[k,s]]/(Log[3^(-k)])},{s,0,500},AxesLabel->{"alpha","f(alpha)"}]

Quit[ ];

Chapter 16: Chaos Control and Synchronization

Chaos Control in the Logistic Map (Figure 16.3(b))

mu=4;switchon=60;itermax=200;x[0]=0.6;k=0.217;
f[x_]=mu x (1-x);

Do[{x[i+1]=f[x[i]]},{i,0,switchon-1}]
Do[{x[i+1]=k f[x[i]],x[i+2]=f[x[i+1]]},{i,switchon,itermax,2}]
p1=ListPlot[Table[{i,x[i]},{i,0,itermax}],PlotJoined->True,PlotRange->All];
p2=ListPlot[Table[{i,x[i]},{i,0,itermax}],PlotStyle->PointSize[0.015]];
Show[{p1,p2}]

Chaos Control in the Hénon Map (Figure 16.6(a))

alpha=1.2;beta=0.4;
k1=-1.8;k2=1.2;

xstar=0.8358;ystar=xstar;

itermax=400;x[0]=0.5;y[0]=0.6;r[0]=(x[0])^2+(y[0])^2;switch=196;
Henon[{x_,y_}]={alpha+beta y-x^2,x};
ControlHenon[{x_,y_}]={(-k1 (x-xstar)-k2 (y-ystar)+alpha)+beta y-x^2,x};

Do[{x[i+1],y[i+1]}=Henon[{x[i],y[i]}],{i,0,switch-1}]
Do[{x[i+1],y[i+1]}=ControlHenon[{x[i],y[i]}],{i,switch,itermax}]
Do[r[i+1]=(x[i+1])^2+(y[i+1])^2,{i,0,itermax}]
p1=ListPlot[Table[{i,r[i]},{i,0,itermax}],PlotJoined->True,PlotRange->All];
p2=ListPlot[Table[{i,r[i]},{i,0,itermax}],PlotStyle->PointSize[0.015]];
Show[{p1,p2},PlotRange->All]

Chaos Synchronization Between Two Lorenz Systems (Figure 16.7(b))

sigma=16;b=4;r=45.92;

Lorenz=NDSolve[{x1'[t]==sigma (x2[t]-x1[t]),x2'[t]==-x1[t] x3[t]+r x1[t]-x2[t],x3'[t]==x1[t] x2[t]-b x3[t],y2'[t]==-x1[t] y3[t]+r x1[t]-y2[t],y3'[t]==x1[t] y2[t]-b y3[t],x1[0]==15,x2[0]==20,x3[0]==30,y2[0]==10,y3[0]==20},{x1,x2,x3,y2,y3},{t,0,100},MaxSteps->Infinity];

ParametricPlot3D[Evaluate[{x1[t],x2[t],x3[t]}/.Lorenz],{t,0,100},PlotPoints->10000]

ParametricPlot[Evaluate[{x3[t],y3[t]}/.Lorenz],{t,50,100},PlotPoints->10000,AxesLabel->{"(Subscript[x, 3])(t)","(Subscript[y, 3])(t)"}]

Generalized Synchronization: The Auxiliary Approach (Figure 16.8(a))

mu=5.7;sigma=16;b=4;r=45.92;g=4; (*No Synchronization*)

RosslerLorenz=NDSolve[{x1'[t]==-(x2[t]+x3[t]),x2'[t]==x1[t]+0.2 x2[t],x3'[t]==0.2+x3[t] (x1[t]-mu),y1'[t]==sigma (y2[t]-y1[t])-g(y1[t]-x1[t]),y2'[t]==-y1[t] y3[t]+r y1[t]-y2[t],y3'[t]==y1[t] y2[t]-b y3[t],z1'[t]==sigma (z2[t]-z1[t])-g(z1[t]-x1[t]),z2'[t]==-z1[t] z3[t]+r z1[t]-z2[t],z3'[t]==z1[t] z2[t]-b z3[t],x1[0]==2,x2[0]==-10,x3[0]==44,y1[0]==30,y2[0]==10,y3[0]==20,z1[0]==31,z2[0]==11,z3[0]==22},{x1,x2,x3,y1,y2,y3,z1,z2,z3},{t,0,200},MaxSteps->Infinity];

ParametricPlot[Evaluate[{y2[t],z2[t]}/.RosslerLorenz],{t,50,200},PlotPoints->10000,AxesLabel->{"(Subscript[y, 2])(t)","(Subscript[z, 2])(t)"}]

mu=5.7;sigma=16;b=4;r=45.92;g=8; (*Synchronization*)

RosslerLorenz=NDSolve[{x1'[t]==-(x2[t]+x3[t]),x2'[t]==x1[t]+0.2 x2[t],x3'[t]==0.2+x3[t] (x1[t]-mu),y1'[t]==sigma (y2[t]-y1[t])-g(y1[t]-x1[t]),y2'[t]==-y1[t] y3[t]+r y1[t]-y2[t],y3'[t]==y1[t] y2[t]-b y3[t],z1'[t]==sigma (z2[t]-z1[t])-g(z1[t]-x1[t]),z2'[t]==-z1[t] z3[t]+r z1[t]-z2[t],z3'[t]==z1[t] z2[t]-b z3[t],x1[0]==2,x2[0]==-10,x3[0]==44,y1[0]==30,y2[0]==10,y3[0]==20,z1[0]==31,z2[0]==11,z3[0]==22},{x1,x2,x3,y1,y2,y3,z1,z2,z3},{t,0,500},MaxSteps->Infinity];

ParametricPlot[Evaluate[{y2[t],z2[t]}/.RosslerLorenz],{t,50,500},PlotPoints->10000,AxesLabel->{"(Subscript[y, 2])(t)","(Subscript[z, 2])(t)"}]

Quit[ ];

Chapter 17: Neural Networks

Generalized Delta Learning Rule (Figure 17.7)

N.B. ENSURE THAT housing.txt IS IN THE CORRECT DIRECTORY. EXPORT A FILE TO SEE WHERE IT IS EXPORTED TO!

X=Import["housing.txt","Table"];
t=X[[All,14]];
tmean=(Max[t]+Min[t])/2;
tstd=(Max[t]-Min[t])/2;
t=N[(t-tmean)/tstd];
Y={X[[All,6]],X[[All,9]],X[[All,13]]};
Y=Transpose[Y];
MnX=List[Mean[Y]];
Std=List[StandardDeviation[Y]];
Dim=Part[Dimensions[Y],1];
ones=Table[{1},{i,Dim}];
X=N[(Y-ones.MnX)/ones.Std];
X=Transpose[X];
X=Prepend[X,Flatten[Table[1,{i,Dim}]]];
SeedRandom[123456];
ww[1]=0.1 Table[{Random[]},{4}];
Epochs=10;
eta=0.001;k=1;
For[n=1,n≤Epochs,n++,For[j=1,j≤Dim,j++,yk=Flatten[Tanh[List[X[[All,j]]].ww[k]]];err=yk-t[[j]];g=Transpose[List[X[[All,j]]]].((1-yk.yk) err);ww[k+1]=ww[k]-eta g;k=k+1]];
For[i=1,i≤4,i++,pts[i]=Table[Flatten[ww[j]][[i]],{j,1,k}]]
For[i=1,i≤4,i++,p[i]=ListPlot[pts[i],PlotStyle->PointSize[0.001]]]
Show[{p[1],p[2],p[3],p[4]},AxesLabel->{"Iterations","Weights"},PlotRange->All]

Example 5 The Discrete Hofield Network

X={{1,1,1,1,1},{1,-1,-1,1,-1},{-1,1,-1,1,1}};
W=Transpose[X].X/5-3 IdentityMatrix[5]/5
{{0,-(1/5),(1/5),(1/5),-(1/5)},{-(1/5),0,(1/5),(1/5),(3/5)},{(1/5),(1/5),0,-(1/5),(1/5)},{(1/5),(1/5),-(1/5),0,(1/5)},{-(1/5),(3/5),(1/5),(1/5),0}}
RandomRelist[x_List]:=Block[{n=x,p},Do[p=Random[Integer,{1,i}];n[[{p,i}]]=n[[{i,p}]],{i,Length[x]}];n];
n=RandomRelist[{1,2,3,4,5}]

xinput={-1,-1,1,1,1};
xtest=xinput;
hsgn[{v_,x_}]=Which[v>0,1,v==0,x,v<0,-1];
For[j=1,j≤5,j++,xtest[[n[[j]]]]=hsgn[{W[[n[[j]]]].xtest,xtest[[n[[j]]]]}]];
If[xtest==X[[1]],Print["Net has converged to X1"],If[xtest==X[[2]],Print["Net has converged to X2"],If[xtest==X[[3]],Print["Net has converged to X3"],Print["Iterate again: May have converged to a spurious steady state"]]]]

Quit[ ];

The Minimal Chaotic Neuromodule (Figure 17.8)

b1=-2;b2=3;w11=-20;w21=-6;w12=6;

Itermax=10000;

Neuromodule[{x_,y_}]={b1+w11/(1+Exp[-x])+w12/(1+Exp[-y]),b2+w21/(1+Exp[-x])}

NeuroMap=N[Drop[NestList[Neuromodule,{1,0.2},Itermax],50]];

ListPlot[NeuroMap,AxesLabel->{"x","y"},PlotStyle->{PointSize[0.001]},Axes->None,Frame->True]

Quit[ ];

Bifurcation Diagram for a Bistable Neuromodule (Figure 17.15)

(*Bifurcation diagram.*)
Start=-5;
Maximum=10;b2=3;w11=7;w12=-4;w21=5;
HalfN=4999;phi=0;N1=1+HalfN;Itermax=2 HalfN+1;x[1]=-10;y[1]=-3;

For[n=1,n≤HalfN,n++,{b1=(Start+n Maximum/HalfN),x[n+1]=N[b1+w11 Tanh[x[n]]+w12 Tanh[0.3 y[n]]],y[n+1]=N[b2+w21 Tanh[x[n]]]}];
p1=ListPlot[Table[Flatten[{Start+n*Maximum/N1,x[n]}],{n,1,HalfN}],PlotStyle->{PointSize[0.001],RGBColor[0,1,0]}];
For[n=N1,n≤Itermax,n++,{b1=(Start+2 Maximum-n Maximum/HalfN),x[n+1]=N[b1+w11 Tanh[x[n]]+w12 Tanh[0.3 y[n]]],y[n+1]=N[b2+w21 Tanh[x[n]]]}];
p2=ListPlot[Table[Flatten[{Start+n*Maximum/N1,x[Itermax+1-n]}],{n,1,HalfN}],PlotStyle->{PointSize[0.001],RGBColor[0,1,0]}];
Show[{p1,p2},PlotRange->All,AxesLabel->{"(Subscript[b, 1])","(Subscript[x, n])"},Frame->True,TextStyle->{FontSize->15},AspectRatio->Automatic]

Quit[ ];