MATLAB Code for funtion lyapunov
function [x,Lce,ErrorQ]=lyapunov(F1,F2,Ic,Ntot,LceFlag,Nskip,OrthFlag);
% lyapunov Computes the trajectory and the Lyapunov characteristic
% exponents for discrete systems via the Householder QR Based Method.
% For more information, see "An efficient QR based method for the
% computation of Lyapunov Exponents", by Hubertus F. von Bremen,
% Firdaus E. Udwadia and Wlodek Proskurowski, Physica D, 1997,
to appear.
% Program written by H. von Bremen 10/10/96
%
% [x,Lce,ErrorQ]=lyapunov('xfile','tanmfile',Ic,Ntot,LceFlag,Nskip,OrthFlag)
% computes the trajectory of a discrete dynamical system
which is
% specified in the M-file xfile.m . It also computes
the Lyapunov
% characteristic exponents (Lce's), provided the tangent
maps are
% specified in the M-file tanmfile.m. The trajectory
is computed for
% up to Ntot iterations. The LCE's are computed starting
at the Nskip+1
% iteration of the trajectory, and are computed for Ntot-Nskip
iterations.
%
%
% INPUT:
% xfile String containing
the user-supplied function file that defines
%
the trajectory.
%
CALL: xmap=fun(x,iter) where fun = 'xfile'
%
x - state vector
%
iter - iteration number
%
xmap - returned next trajectory point
%
{x(iter+1)=xmap=fun(x(iter),iter)}
% tanmfile String containing the user-supplied
file that defines the
%
tangent map to the trajectory.
%
CALL: tanmap=fun(x,iter) where fun = 'tanmfile'
%
x - state vector
%
iter - iteration number
%
tanmap - returned next tangent map
%
{Tangent_Map(iter+1)=tanmap=fun(x(iter),iter)}
% Ic
Initial conditions for the trajectory, given as a column
%
vector
% Ntot Number
of iterations (time steps) to be used to calculate
%
the trajectory of the system.
% LceFlag: Flag used to determine if
Lce's are computed.
%
-If LceFlag = 1, then Lce's are computed.
%
-If LceFlag = 0, then LCE's are not computed and only
%
the trajectory is computed.
% Nskip: Number of iterations
to be skipped before the LCE's are
%
to be computed. First Nskip points along the trajectory are
%
skipped in order to allow the disturbances due to the
%
initial conditions to diminish, before LCE's are computed.
% OrthFlag: Flag used to determine if error
in orthogonality is to
%
be computed, and if computed, then of what type.
%
-If OrthFlag = 0 (recommended value), then the error in
%
orthogonality is not computed.
%
-If OrthFlag = 1, then the determinant error in orthogonality
%
is computed. Error in Orthogonality = abs(1-abs(det(Q)))
%
-If OrthFlag = 2, then the Two-Norm error in orthogonality
%
is computed. Error in Orthogonality = norm(Q'*Q-I)
%
NOTE: when using OrthFlag = 1 or 2, the computational cost is
%
increased significantly (using 2 is more expensive than 1).
%
% OUTPUT:
% x
Trajectory of dynamical system given in matrix form, where the
%
k-th column are the trajectory coordinates for the k-th
%
iteration (the first column is Ic). Size(x)=[n,Ntot+1].
% Lce Computed
Lyapunov exponents given in matrix form, where the
%
k-th column are the LCE's for the k-th iteration of LCE's.
%
Size(Lce)=[n,Ntot-Nskip].
% ErrorQ Error in othogonality
given in vector form, where the k-th
%
entry is the error in orthogonality of the k-th Lce iteration.
%
Size(ErrorQ)=[1,Ntot-Nskip].
%
-If OrthFlag = 1, ErrorQ = abs(1-abs(det(Q)))
%
-If OrthFlag = 2, ErrorQ = Norm(Q'*Q-I)
%
%
% **SAMPLE:
% *************
% Sample INPUT:
% Ic=[0.2 0.3]'; Ntot= 1000; LceFlag=1; Nskip=100;
OrthFlag=0;
%
% Sample function files for the trajectory and the tangent
map:
% ---------------------------
% Sample for xfile='map'; coupled logistic maps (stored
in m-file map.m)
% function xmap=map(x,iter)
% d=.2; r=3.6;
% xm(1)=d*r*x(1)*(1-x(1)) +
(1-d)*r*x(2)*(1-x(2));
% xm(2)=(1-d)*r*x(1)*(1-x(1)) + d*r*x(2)*(1-x(2));
% xmap=[xm(1) xm(2)]';
% ---------------------------
% Sample for tanmfile='tanmap'; Jacobian for coupled logistic
maps (stored in
% m-file tanmap.m)
% function [T]=tanmap(x)
% d=.2; r=3.6;
% T= [ d*r*(1-2*x(1))
(1-d)*r*(1-2*x(2))
% (1-d)*r*(1-2*x(1))
d*r*(1-2*x(2)) ];
% ---------------------------
% Note: the above two sample function files are provided
in the files map.m
% and tanmap.m.
% --------------------------------
% *************
% Sample CALL
% [x,Lce,ErrorQ]=lyapunov('map','tanmap',Ic,Ntot,LceFlag,Nskip,OrthFlag);
% *************
% Sample OUTPUT:
% To plot the LCE's for each iteration you can use:
% subplot(2,1,1), plot(Lce(1,:))
% subplot(2,1,2), plot(Lce(2,:))
% To plot the trajectory coordinates for each iteration
you can use:
% subplot(2,1,1), plot(x(1,:))
% subplot(2,1,2), plot(x(2,:))
%
% Note: By setting OrthFlag=1 (or = 2), you can plot the
error in
% orthogonality by using: plot(ErrorQ)
%
% In general we have:
% To plot the k-th LCE for each iteration: plot(Lce(k,:))
% To plot the k-th coordinate of the trajectory:
plot(x(k,:))
% If OrthFlag is 1 or 2,then to plot the error in orthogonality:
plot(ErrorQ)
%
% INPUT ERROR HANDLING
if (LceFlag~=0)&(LceFlag~=1)
disp('You have entered an invalid value for LceFlag')
disp('You should choose LceFlag to be either 0 or 1 (see help
file)')
disp('The value of LceFlag was set to 0 ')
disp(' ')
LceFlag=0;
end
if (OrthFlag~=0)&(OrthFlag~=1)&(OrthFlag~=2)
disp('You have entered an invalid value for OrthFlag')
disp('You should choose OrthFlag to be either 0 or 1 or 2 (see
help file)')
disp('The value of OrthFlag was set to 0')
disp(' ')
OrthFlag=0;
end
if Ntot<Nskip
disp('You have entered an invalid value for either Ntot
or Nskip')
disp('You must have Ntot>=Nskip, (see help file)')
disp('The values of Nskip and Ntot have been switched')
disp(' ')
temp1=Nskip;
Nskip=Ntot;
Ntot=temp1;
end
% INITIALIZATION PART
m=Ntot-Nskip;
n=max(size(Ic));
x=zeros(n,Ntot+1);
% END INITIALIZATION
x(:,1)=Ic;
for j=2:Nskip+1
Ic=feval(F1,Ic);
x(:,j)=Ic;
end % (j)
if LceFlag==0
for j=1:m
Ic=feval(F1,Ic);
x(:,j+Nskip+1)=Ic;
end % (j)
else
r=zeros(1,n);gama=zeros(1,n);
q=eye(n);sl=zeros(n,1);Lce=zeros(n,m);
if OrthFlag==1|OrthFlag==2
ErrorQ=zeros(1,m);
end % (if)
A=feval(F2,Ic);
if OrthFlag==0
a=A;
end;
N2=Nskip+1;
for i=1:m
if OrthFlag~=0
a=A*q;
Ic=feval(F1,Ic);
x(:,i+N2)=Ic;
A=feval(F2,Ic);
B=eye(n);
elseif OrthFlag == 0
Ic=feval(F1,Ic);
x(:,i+N2)=Ic;
B=feval(F2,Ic);
end % (if OrthFlag)
% Computation of the Factorization
for k=1:n-1
% computation of the reflector
if a(k,k)<0
b=-1;
else
b=1;
end
sig=sqrt(a(k:n,k)'*a(k:n,k));
gama=sig*(sig+abs(a(k,k)));
r(k)=-b*sig;
a(k,k)=a(k,k)-r(k);
% end computation of the reflector
p=k+1;
for j=p:n
b=a(k:n,k)'*a(k:n,j)/gama;
a(p:n,j)=a(p:n,j)-b*a(p:n,k);
end
% Computation of B*Q
for j=1:n
b=B(j,k:n)*a(k:n,k)/gama;
B(j,k:n)=B(j,k:n)-a(k:n,k)'*b;
end
end
r(n)=a(n,n);
% End Computation of the factorization
if OrthFlag~=0
q=B;
if OrthFlag==1
ErrorQ(i)=abs(1-abs(det(q)));
elseif OrthFlag==2
ErrorQ(i)=norm(q'*q-eye(n));
end % (if)
elseif OrthFlag == 0 % (if OrthFlag)
a=B;
end % (if OrthFlag)
sl=sl+log(abs(r'));
Lce(:,i)=sl/i;
end % (i)
end % (else)
The University of Southern California does not screen or control the content on this website and thus does not guarantee the accuracy, integrity, or quality of such content. All content on this website is provided by and is the sole responsibility of the person from which such content originated, and such content does not necessarily reflect the opinions of the University administration or the Board of Trustees