Using MATLAB to visualize Chirkov Standard Map

Question

I am looking to visualize the Chirikov Standard Map: $$ \theta_{n+1} = \left( \theta_n + I_n - \frac{K}{2\pi} \sin(2\pi\theta_n)\right) \mod 1$$ $$ I_{n+1} = I_n -\frac{K}{2\pi}\sin(2\pi\theta_n)$$ where $K=1$ using MATLAB. I have the following code, but the plot generated does not seem right given my understanding of the chaotic behavior of the Chirikov Standard Map. My initial condition is the unit square.

close all; 
clear all; 
clc; 

% Initial Condition: 
% Top: 
xTop = linspace(0,1,10); 
xTop = xTop'; 
yTop = ones(length(xTop),1)*1; 
top = []; 
top = [top xTop];
top = [top yTop]; 
% Bottom: 
xBottom = linspace(0,1,10); 
xBottom = xBottom'; 
yBottom = ones(length(xBottom),1)*0; 
bottom = []; 
bottom = [bottom xBottom];
bottom = [bottom yBottom]; 
% Right: 
yRight = linspace(0,1,10); 
yRight = yRight'; 
xRight = ones(length(yRight),1)*1; 
right = []; 
right = [right xRight]; 
right = [right yRight]; 
% Left: 
yLeft = linspace(0,1,10); 
yLeft = yLeft'; 
xLeft = ones(length(yLeft),1)*(0); 
left = []; 
left = [left xLeft]; 
left = [left yLeft]; 

thetaCurr = []; 
thetaCurr = [thetaCurr xTop xBottom xRight xLeft]; 
iCurr = []; 
iCurr = [iCurr yTop yBottom yRight yLeft];
iterations = 2; 
hold on; 

% Forward Iteration of the Chirkov Standard Map: 
for t=1:iterations
   thetaNext = thetaCurr+iCurr-(1/(2*pi))*sin(2*pi*thetaCurr); 
   thetaNext = mod(thetaNext,1); 
   iNext = iCurr-(1/(2*pi))*sin(2*pi*thetaCurr); 
   for i=1:4 
      p2 = plot(thetaNext(:,i),iNext(:,i),'r');  
   end
   iCurr = iNext; 
   thetaCurr = thetaNext; 
end

Any and all help is greatly appreciated!

Answer 1

Probably the problem may be with the value of constant 'K'.It should be above 18 for Chirikov map to behave chaotically.

Answer 2

The iteration code is working quite well, and the plot can be improved by plotting marks instead of lines, and increasing the number of iterations to some extent. Here is a code with K=0.8, where I kept the same notation but simplified the initial condition by using the "meshgrid" function. Hope this is useful (although with a few years of delay..)

clc; 

x=0:0.1111:1;
y=0:0.1111:1;
[thetaCurr,iCurr]=meshgrid(x,y);

iterations = 100;
K= 0.8  % 0.971635

hold on; 

for t=1:iterations
   thetaNext = thetaCurr+iCurr-(K/(2*pi))*sin(2*pi*thetaCurr); 
   thetaNext = mod(thetaNext,1); 
   iNext = iCurr-(K/(2*pi))*sin(2*pi*thetaCurr); 
   for i=1:size(x,2) 
      p2 = plot(thetaNext(:,i),iNext(:,i),'r*',"MarkerSize",1);  
   end
   iCurr = iNext; 
   thetaCurr = thetaNext; 
end