%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Example Octave/MATLAB programs to produce the graphs %% for Math 3520, HW 6, Spring 2026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Ex.6 %% Henon area-preserving mapping N = 20000; % Number of points to plot %% Various alphas %%% First alpha, four initial points % Henon's choice cosalpha = 0.24; sinalpha = sqrt(1-cosalpha*cosalpha); % parameters henonapm = @(X) [X(1)*cosalpha - (X(2)-X(1)*X(1))*sinalpha, X(1)*sinalpha + (X(2)-X(1)*X(1))*cosalpha]; % see p.488 X=[0.5,0.5]; % Nonzero initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=0.24, X0=[0.5,0.5]"); xlabel("Xn"); ylabel("Yn"); X(1,:)=[0.57,0.16]; % Suggested initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=0.24, X0=[0.57,0.16]"); xlabel("Xn"); ylabel("Yn"); X(1,:)=[0.7,0.7]; % Near-edge Initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=0.24, X0=[0.7,0.7]"); xlabel("Xn"); ylabel("Yn"); X(1,:)=[-0.4,0.1]; % Other Initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=0.24, X0=[-0.4,0.1]"); xlabel("Xn"); ylabel("Yn"); %%% Second alpha, same four initial points % Larger cos(alpha) cosalpha = 0.67; sinalpha = sqrt(1-cosalpha*cosalpha); % parameters henonapm = @(X) [X(1)*cosalpha - (X(2)-X(1)*X(1))*sinalpha, X(1)*sinalpha + (X(2)-X(1)*X(1))*cosalpha]; % see p.488 X=[0.5,0.5]; % Nonzero initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=0.67, X0=[0.5,0.5]"); xlabel("Xn"); ylabel("Yn"); X(1,:)=[0.57,0.16]; % Suggested initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=0.67, X0=[0.57,0.16]"); xlabel("Xn"); ylabel("Yn"); X(1,:)=[0.7,0.7]; % Near-edge Initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=0.67, X0=[0.7,0.7]"); xlabel("Xn"); ylabel("Yn"); X(1,:)=[-0.4,0.1]; % Other Initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=0.67, X0=[-0.4,0.1]"); xlabel("Xn"); ylabel("Yn"); %%% Third alpha, same four initial points % Smaller negative cos(alpha) cosalpha = -0.17; sinalpha = sqrt(1-cosalpha*cosalpha); % parameters henonapm = @(X) [X(1)*cosalpha - (X(2)-X(1)*X(1))*sinalpha, X(1)*sinalpha + (X(2)-X(1)*X(1))*cosalpha]; % see p.488 X=[0.5,0.5]; % Nonzero initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=-0.17, X0=[0.5,0.5]"); xlabel("Xn"); ylabel("Yn"); X(1,:)=[0.57,0.16]; % Suggested initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=-0.17, X0=[0.57,0.16]"); xlabel("Xn"); ylabel("Yn"); X(1,:)=[0.7,0.7]; % Near-edge Initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=-0.17, X0=[0.7,0.7]"); xlabel("Xn"); ylabel("Yn"); X(1,:)=[-0.4,0.1]; % Other Initial values [x0,y0] for(i=2:N) X(i,:) = henonapm(X(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([-1,1]); ylim([-1,1]); % square [-1,1]x[-1,1] plot(X(:,1), X(:,2), "."); % plot trajectory title("Henon Area-Preserving Map, cos(alpha)=-0.17, X0=[-0.4,0.1]"); xlabel("Xn"); ylabel("Yn"); %%% Ex. 7b %% Standard mapping template: % ynp1=mod(yn+k*sin(xn), 2*pi); % compute y_{n+1} first % xnp1=mod(xn+ynp1, 2*pi); % ...then compute x_{n+1} % ... or use substitution in the anonymous function % stdmap = @(X)[mod(X(1)+X(2)+ k*sin(X(1),2*pi), mod(X(2)+k*sin(X(1)),2*pi)]; %% Initialization valid for all graphs N = 20000; % Number of points to plot twopi = 2*pi; % Precompute 2*pi for efficiency %% Various ks %%% First k, four initial points k=0; stdmap = @(X)[mod(X(1)+X(2),twopi), mod(X(2),twopi) ]; Xa=[1,1]; % first initial point for(i=2:N) Xa(i,:) = stdmap(Xa(i-1,:)); end Xb=[1,5]; % second initial point for(i=2:N) Xb(i,:) = stdmap(Xb(i-1,:)); end Xc=[2,twopi/17]; % third initial point for(i=2:N) Xc(i,:) = stdmap(Xc(i-1,:)); end Xd=[twopi/3,pi]; % fourth initial point for(i=2:N) Xd(i,:) = stdmap(Xd(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([0,twopi]); ylim([0,twopi]); % square S plot(Xa(:,1), Xa(:,2), "."); % first orbit plot(Xb(:,1), Xb(:,2), "."); % second orbit plot(Xc(:,1), Xc(:,2), "o"); % third orbit plot(Xd(:,1), Xd(:,2), "x"); % fourth orbit title("Standard Map, k=0, X0=[1,1], [1,5], [2,2pi/17], [2pi/3,pi]"); xlabel("Xn"); ylabel("Yn"); %%% Second k, same four initial points k=0.5; stdmap = @(X)[mod(X(1)+X(2)+k*sin(X(1)),twopi), mod(X(2)+k*sin(X(1)),twopi) ]; Xa=[1,1]; % first initial point for(i=2:N) Xa(i,:) = stdmap(Xa(i-1,:)); end Xb=[1,5]; % second initial point for(i=2:N) Xb(i,:) = stdmap(Xb(i-1,:)); end Xc=[2,twopi/17]; % third initial point for(i=2:N) Xc(i,:) = stdmap(Xc(i-1,:)); end Xd=[twopi/3,pi]; % fourth initial point for(i=2:N) Xd(i,:) = stdmap(Xd(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([0,twopi]); ylim([0,twopi]); % square S plot(Xa(:,1), Xa(:,2), "."); % first orbit plot(Xb(:,1), Xb(:,2), "."); % second orbit plot(Xc(:,1), Xc(:,2), "."); % third orbit plot(Xd(:,1), Xd(:,2), "."); % fourth orbit title("Standard Map, k=0.5, X0=[1,1], [1,5], [2,2pi/17], [2pi/3,pi]"); xlabel("Xn"); ylabel("Yn"); %%% Third k, same four initial points k=1; stdmap = @(X)[mod(X(1)+X(2)+k*sin(X(1)),twopi), mod(X(2)+k*sin(X(1)),twopi) ]; Xa=[1,1]; % first initial point for(i=2:N) Xa(i,:) = stdmap(Xa(i-1,:)); end Xb=[1,5]; % second initial point for(i=2:N) Xb(i,:) = stdmap(Xb(i-1,:)); end Xc=[2,twopi/17]; % third initial point for(i=2:N) Xc(i,:) = stdmap(Xc(i-1,:)); end Xd=[twopi/3,pi]; % fourth initial point for(i=2:N) Xd(i,:) = stdmap(Xd(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([0,twopi]); ylim([0,twopi]); % square S plot(Xa(:,1), Xa(:,2), "."); % first orbit plot(Xb(:,1), Xb(:,2), "."); % second orbit plot(Xc(:,1), Xc(:,2), "."); % third orbit plot(Xd(:,1), Xd(:,2), "."); % fourth orbit title("Standard Map, k=1, X0=[1,1], [1,5], [2,2pi/17], [2pi/3,pi]"); xlabel("Xn"); ylabel("Yn"); %%% Fourth k, same four initial points k=2; stdmap = @(X)[mod(X(1)+X(2)+k*sin(X(1)),twopi), mod(X(2)+k*sin(X(1)),twopi) ]; Xa=[1,1]; % first initial point for(i=2:N) Xa(i,:) = stdmap(Xa(i-1,:)); end Xb=[1,5]; % second initial point for(i=2:N) Xb(i,:) = stdmap(Xb(i-1,:)); end Xc=[2,twopi/17]; % third initial point for(i=2:N) Xc(i,:) = stdmap(Xc(i-1,:)); end Xd=[twopi/3,pi]; % fourth initial point for(i=2:N) Xd(i,:) = stdmap(Xd(i-1,:)); end figure; hold on; % New graph, accumulate points xlim([0,twopi]); ylim([0,twopi]); % square S plot(Xa(:,1), Xa(:,2), "."); % first orbit plot(Xb(:,1), Xb(:,2), "."); % second orbit plot(Xc(:,1), Xc(:,2), "."); % third orbit plot(Xd(:,1), Xd(:,2), "."); % fourth orbit title("Standard Map, k=2, X0=[1,1], [1,5], [2,2pi/17], [2pi/3,pi]"); xlabel("Xn"); ylabel("Yn"); %%% Ex.9 % Henon map trapping region % Specified a,b for this exercise: a=1.4; b=0.3; % Henon mapping formula, works with N rows of 2-vectors henon = @(X)[X(:,2)+1-a*X(:,1).^2, b*X(:,1)]; % Parametric sides of the quadrilateral A=[-1.33,0.42]; B=[1.32,0.133]; C=[1.245,-0.14]; D=[-1.06,-0.5]; ts = 0:0.01:1; % 100 steps of the t parameter % (a) Plot the 4 sides AB, BC, CD, DA figure; hold on; AB=ts'*(B-A)+A; BC=ts'*(C-B)+B; CD=ts'*(D-C)+C; DA=ts'*(A-D)+D; QXY = [AB;BC;CD;DA]; TQ = henon(QXY); plot(QXY(:,1),QXY(:,2)); plot(TQ(:,1),TQ(:,2)) title("Ex.9a: Quadrilateral ABCD (Blue) and T(ABCD) (Red)"); xlabel("Xn"); ylabel("Yn"); % (b) Direction vectors and normals vAB = B-A vBC = C-B vCD = D-C vDA = A-D %%%%% Ex.10c %%%% Runge-Kutta 4th order solver in 3D: function X = rk4b(X0,F,steps,h) % Input % X0 = [x0,y0,z0] Row vector, initial position % F = function handle, computes Xdot(x,y,z) % steps = number of time steps to compute (e.g., 100) % h = time increment per step (e.g., 1/steps ) % Output % X = (1+steps)x(3) matrix of (x,y,z) values along % the trajectory starting at X0 X = zeros(1+steps,3); % Solution to be computed X(1,:) = X0; % Initial point for n=1:steps f1 = F(X(n,:)); f2 = F(X(n,:)+f1*h/2); f3 = F(X(n,:)+f2*h/2); f4 = F(X(n,:)+f3*h); X(n+1,:) = X(n,:) +h*(f1+2*f2+2*f3+f4)/6; end end %%%% From Strogatz, 3rd edition, p.472 X0=[1,1,1]; % same as Lorenz's initial value a=0.2; % from Strogatz b=0.2; % from Strogatz, same b=a c=5.7; % from Strogatz, to get chaos Xdot = @(V) [-V(2)-V(3), V(1)+a*V(2), b+V(3).*(V(1)-c)]; %%%% Trajectories h=0.01; % Use same time step for same accuracy X = rk4b(X0, Xdot, 20000,h); % Complete trajectory to t=200. figure; hold on; plot3(X(:,1),X(:,2),X(:,3)); % initial search for x,y,z,limits xlabel("X"); ylabel("Y"); zlabel("Z"); title("Trajectory for the Rossler System") % Conjecture: trapped in -15