Back to Gallery
  • Remix
  • Share
  • New Entry

  • leon t

  • /

  • Stable Manifold of Lorenz System

on 28 Nov 2023
  • 22
  • 134
  • 0
  • 2
  • 1915
drawframe(1);
Write your drawframe function below
function drawframe(f)
b=8/3;r=28;s=10;
L=@(t,v)[s*(v(2)-v(1));
v(1).*(r-v(3))-v(2);
v(1).*v(2)-b.*v(3)];
Vss=[-((2209-9*1201^(1/2))^(1/2)*(9*122^(1/2)+146522^(1/2)))/34160
(122^(1/2)*(2209-9*1201^(1/2))^(1/2))/610
0];
Vs=[0
0
1];
AL=150;
r0=1;
NA=10;
AcL=linspace(0,AL,NA+1);AcL(1)=[];
R=r0+AL;
N=801;
TS=0:-1e-3:-30;
HP=sym(pi)/2;
q=444e-18;
Ag=[sym([0
223e-23
668e-23
178e-22
223e-21
89e-19
89e-18
371e-16
542e-15
153e-13
582e-12
102e-10])-HP
-1.570796
-1.570766
-1.570324
-0.237841
1.556471
1.570748
1.570794
HP-sym(1e-10*[1790
549
156
66
31.4
22.1
11.3
4.86
2.57
1.95
1.66
1.46])
1.57079632665462698426+cumsum(sym([0
2.080558e-13
509e-16
182e-16
311e-17
3*q/2
q
q
q
q
q
q
3*q/2
2*q
q
q
q
q
q
2*q
289e-17
16e-15]))
HP-sym(1e-10*[1.37
1.28
1.16
1.02
.902
.795
.746
.689
.632
.574
.517
.458
.403
.353
.294
.235
.199
.174
.146
.117
.0874
.053
0])];
Lc=nan(numel(Ag),(3*NA+3));
for i=1:numel(Ag)
CP=double(r0*(Vss*cos(Ag(i))+Vs*sin(Ag(i))));
Lc(i,:)=[reshape(CP,[1,3]),LocusCalculation(CP,L,AcL,R,TS,N)];
end
j=3*NA+1:-3:1;
X=[Lc(:,[j,1]);...
-flipud(Lc(:,[j,1]))];
Y=[Lc(:,[j+1,2]);...
-flipud(Lc(:,[j+1,2]))];
Z=[Lc(:,[j+2,3]);...
flipud(Lc(:,[j+2,3]))];
[u,v]=size(X);
mC=ones(u,1)*linspace(0,AL,v);
j=1:min(f*4+1,u);
hold on
j=1:min(f*4+1,u);
surf(X(j,:),Y(j,:),Z(j,:),mC(j,:),'FaceAlpha',3/4,'MeshStyle','row','LineWidth',.2)
colormap(flipud(hot(u)))
plot3(X(j,1),Y(j,1),Z(j,1),'b-','LineWidth',2)
xlabel x
ylabel y
zlabel z
title('Stable Manifold of Lorenz System')
axis(100*[-1 1 -1 1 -2 5/2])
axis off
daspect(ones(1,3))
view([-152+9*f,0])
if f>47
text(-30,40,180,sprintf(' TouAkira\n@iLoveMatlab.CN'),'FontWeight','bold');
end
function O=LocusCalculation(SP,L,AcL,R,TS,N)
Sl=ode45(L,TS,SP,odeset('Events',@(t,y)Evt(t,y,R)));
if numel(Sl.x)<N
if abs(Sl.x(end))<abs(TS(end))
V=linspace(Sl.x(1),Sl.x(end),N);
else
V=linspace(TS(1),TS(end),N);
end
else
V=TS;
end
[~,Sp]=deval(Sl,V);
Sq=vecnorm(Sp,2,1);
A=cumtrapz(abs(V),Sq);
T=sign(A-AcL(:));
G=arrayfun(@(i)find(T(i,:)==-1|T(i,:)==0,1,'last'),1:size(T,1))';
GT=(V(G+1)+V(G))/2;
S=deval(Sl,GT);
O=S(:)';
end
function[V,T,D]=Evt(t,y,R)
V=norm(y(:,end))-R;
T=1;
D=0;
end
end
Animation
Remix Tree