% HW3: Henon bifurcation diag for Compu Expt 2.2 p.76
% barnett 10/15/09
N = 300; % # iters (any more and is slow)
x = zeros(2,N+1); b = -0.3;
figure; disp('please wait about 8 seconds...')
for a=0:0.01:2.2
f = @(x) [a-x(1,:).^2+b*x(2,:); x(1,:)]; % Henon map w/ params a,b
x(:,1) = [.5;.5]; % initial value, better than (0,0)
for n=1:N, x(:,n+1) = f(x(:,n)); end % iterate
plot(a, x(1, end-100:end), 'k.'); hold on; % only plot the final 100 iters
end
% Now, for some fun: want to make it faster? Vectorize by evolving all
% the different cases of a at once! We avoid explicit (thus slow) double loop:
clear
tic; % start the timer
N = 300; b = -0.3;
a = 0:0.01:2.2;
x = [.5;.5] * ones(size(a));
figure;
for n=1:N, x = [a-x(1,:).^2+b*x(2,:); x(1,:)]; % note: a is list
if n>=N-100, plot(a, x(1,:), 'k.'); hold on; end % add to plot per iteration
end
toc % elapsed time is now 0.27 secs, a factor 30 faster than other method!
% you may even decrease a step to 0.001 and it still takes only 0.8 sec