% global unstable manifold of the origin in the Lorenz system s = 10; rh = 28; b = 8/3; % Lorenz system dim = 3; v = @(x) [s*(x(:,2)-x(:,1)) ... rh*x(:,1)-x(:,2)-x(:,1).*x(:,3) ... x(:,1).*x(:,2)-b*x(:,3)]; h = 0.01; n = 10; f = @(x) rk4(v,x,h,n); % flow map n = 20; X1 = linspace(-1,1,n)'; E = ones(size(X1)); X = [ X1 -E -E; X1 -E E; X1 E -E; X1 E E; ... -E X1 -E; -E X1 E; E X1 -E; E X1 E; ... -E -E X1; -E E X1; E -E X1; E E X1]; % sample points c = [0 0 0]; r = 2*[25 25 25]; tree = Tree(c,r); % the box collection %% continuation algorithm newBox = 1; todo = 2; % flag values x = [0 0 0]; % equilibrium depth = 21; % refinement level tree.insert(x'*ones(1,size(X,1))+1E-6*X', depth, newBox); % insert boxes touched by a small neighborhhod of the equilibrium n0 = 0; n1 = tree.count(depth); while n1 > n0 % while boxes are being added tree.change_flags('all', newBox, todo, depth); % need to sample from new boxes b = tree.first_box(depth); while (~isempty(b)) flag = b(2*dim+1); if(bitand(flag, todo)) c = b(1:dim); r = b(dim+1:2*dim); % center and radius of current box p = X*diag(r) + ones(size(X))*diag(c); % sample points in current box tree.insert(f(p)', depth, newBox, 0); % map points and insert hit boxes end b = tree.next_box(depth); end tree.unset_flags('all', todo, depth); n0 = n1; n1 = tree.count(depth); fprintf('%d boxes\n', n1); clf; boxplot3(tree,'alpha',0.9); axis tight; light('position',[-1,2,2]); view(20,30); drawnow; end %% cleanup delete(tree);