# This code computes the Pecora-Carroll exponents for a real eigenvalue # You could use a different chaotic system if you wanted # rossler attractor f(x,y,z)=-y-z g(x,y,z)=x+a*y h(x,y,z)=b*x-c*z+x*z par a=.36,b=.4,c=4.5 init x=0,y=-4.3,z=-.054 # these are parameters for the numerical method par eps=.01,dt=.05 # # coupling stuff par mu=0,kx=1,ky=0,kz=0 # evaluate the right-hand sides xdot=f(x,y,z) ydot=g(x,y,z) zdot=h(x,y,z) # perturb along the preferred direction xx=x+eps*vx yy=y+eps*vy zz=z+eps*vz # Note that [F(X + eps V)-F(x)]/eps = A V + O(eps) # so we compute A V without using the actual Jacobian # we also add the coupling: # K = diag(kx,ky,kz) so we get # (A - mu K) V mu ios the overall strength of coupling # we integrate this using Euler dx=vx+dt*((f(xx,yy,zz)-xdot)/eps-mu*kx*vx) dy=vy+dt*((g(xx,yy,zz)-ydot)/eps-mu*ky*vy) dz=vz+dt*((h(xx,yy,zz)-zdot)/eps-mu*kz*vz) # now we figure out the magnitude of the change in length amp=sqrt(dx^2+dy^2+dz^2) # # eulers method for integration of the Rossler equations x'=x+dt*xdot y'=y+dt*ydot z'=z+dt*zdot # update the unit vector (vx,vy,vz) vx'=dx/amp vy'=dy/amp vz'=dz/amp # accumulate the log of the expansion/compression ss'=ss+log(amp) # we average it (note that time is discrete so multiply t by dt) aux liap=ss/max(dt,t*dt) # keep track of mu aux muu=mu # initialize (vx,vy,vz) as a unit vector init vx=1 ##### # integrate 80000*dt = 4000 time steps # only keep the last one (trans=79999) @ total=80000,meth=discrete,nout=50 @ trans=79999 @ maxstor=1000000,bound=100000 # plot set up @ xp=muu,yp=liap,xlo=0,xhi=4,ylo=-2,yhi=1 done