chaos.m

混沌识别

APPENDIX
I. Dynamical Systems Analysis—MATLAB Programs
A Numerical Integration of Three Coupled ODEs
% File: VanDerPolSolv.m
%
% Numerical integration of Duffing-Van Der Pol Oscillator.
% Ordinary differential equation solver using in-built
% Matlab function (ode45).
%
% Uses: VanDerPol.m
% Set initial conditions for the three differential equations
xO=[l;O;O];
% Set integration time
timePoints=5000:0.3:7500;
% Solve equations
options=odeset('RelTol',le-10,'AbsTol',le-12);
[t,x]=ode45('VanDerPol',timePoints,x0);
function xp-VanDerPol(t,x)
%
% Duffing-Van Der Pol equations expressed as a first order system
% (see equations 3-5). Parameter values (mu, f and omega)
% can be chosen to demonstrate various trajectories (see text).
%
% IN:
% t: time (not used, but necessary for ODE solver)
% x: Input vector
%
% OUT:
% xp: dx/dt
% Current parameters demonstrate chaotic behavior of the oscillator.
mu=0.2; f=1.0; omega=0.94;
% define equations
xp(l)=x(2);
xp(2)=mu*(l-x(l)*2)*x(2)-x(l)A3+f*cos(x(3));
xp(3)=omega;
% transpose into column vector for ODE solver
xp=xp';
B. Three-Dimensional Phase Space Plots
% File: PSplot.m
%
% Three-dimensional phase space plot of three coupled ODEs
% (data assumed to be stored in matrix 'x' from VanDerPolSolve.m).
%
% Transform x3 to sin(x3)
xnew=x;
xnew(:,3)=sin(x(:,3)) ;
% Generate 3D plot
plot3(xnew(:,1),xnew(:,2),xnew(:,3));
xlabel('xl') ;
ylabeK fx2' ) ;
zlabeK 'sin(x3) ' )
title('Duffing-Van der Pol oscillator.');
rotate3d on
view( [-5,58]) ;
30 Chapter 1 Nonlinear Dynamics Time Series Analysis
C. Two-Dimensional Surface of Section
% File: SurfaceOfSection.m
%
% Two dimension Poincare surface of section on plane x3=0
% (data assumed to be stored in matrix 'xnew' generated from PSplot.m).
%
% Start with empty surface
clear Surface
OldDistance=0;
k=l;
for i=l:size(xnew)*[l;0]
NewDistance=xnew(i,3);
if (NewDistance>=0 * OldDistance<0)
% Add new point to the surface
TotalDistance=NewDistance-OldDistance;
Surface(k,:)=xnew(i-l,:)-(OldDistance/Total Distance)* ...
(xnew(i,:)-xnew(i-1,:));
k=k+l;
end
OldDistance=NewDistance;
end
% Generate 2D plot
plot(Surface(:,1),Surface(:,2),'*');
xlabeK 'xl') ;
ylabeK 'x2') ;
title('Poincare Surface of Section.');
D. Lyapunov Exponents for Three Coupled ODEs
% File: LyapunovSolver.m
%
% Calculate all the Lyapunov exponents for the Duff ing-Van der Pol oscillator .
%
% Uses: IntegrateVDPSystem.m, GramSchmidt.m
% Step 1: Construct a unit hypersphere (Dimension=Dim; axis at time 0:
% a(l)=[l;0;0...;0], a(2)=[0;1;0...;0],...,
a(n)=[0;0;0...;l]
% centered on an initial point (x(l,:))
% n=number of iterations to average
n=5000;
% tau
tau=0.1;
Sigma=zeros(n,3) ;
% Initial conditions
x0=[l;0;0];
a0=eye(3);
x=zeros(3,n);
% j=l (Iteration variable)
for j=l:n
% Step2: Integrate the nonlinear equation of motion over a
% characteristic time scale tau to obtain x(j*tau).
% Integrate the variational equations to obtain the evolved
% axis vectors (al,a2,...,an) for time tau.
[xpoint,a]=IntegrateVDPSystem(xO,aO,tau);
% x: numerical solution of the system {NumberOfPoints.Dim}
x(:,j)=xpoint;
% Step3: Apply Gram-Schmidt re-orthonormalization to the axis vectors.
[aUnit]=GramSchmidt(a);
%Step4: Compute Sigma where
% Sigma(j,l)=log(norm(a(l)),
Sigma(j,2)=log(norm(a(2)),...,
Appendix 31
% Sigma(j,n)=log(norm(a(n))
for i=l:3
aOrth(:,i)=dot(a(:,i),aUnit(: ,i))*aUnit(:, i ) ;
end
for i=l:3
Sigma(j,i)=log(norm(aOrth(:, i),2));
end
xO=xpoint;
aO=aUnit;
end
% Step 5: Compute the Lyapunov exponents
Lyapunov=(l/(n*tau))*sum(Sigma,1);
function [x, a]=IntegrateVDPSystem(xO,aO,tau)
%
% This function integrates the Duffing-Van der Pol equations
% and the variational equations returning the final result at time tau
% for initial conditions xO and aO (3 vectors associated with variational
% equations).
%
% IN:
% xO: Vector {3} of initial conditions
% aO: Matrix {3,3} with initial vectors associated with variational equations.
% tau: Final time.
%
% OUT:
% x: Vector {3} with solution at time tau
% a: Matrix {3,3} with solution of 3 vectors at time tau.
%
% Uses: VanDerPolVarEqu.m
sO=zeros(12,l);
sO=xO;
for i=l:3
sO(i*3+l:i*3+3)=aO(:,i);
end
to=o
tfinal=tO+tau;
options=odeset('RelTol', le-10,'AbsTol',le-12);
[t,s]=ode45('VanDerPolVarEqu',[tO tfinal], sO);
DataLength=size(s)*[l;0];
x=s(DataLength,1:3)';
for i=l:3
a( : ,i)=s(DataLength,i*3+l:i*3+3) ' ;
end
function sp=VanDerPolVarEqu(t,s)
%
% Duffing-Van Der Pol equations expressed as a first order system
% (see equations 3-5). Parameter values (mu, f and omega)
% can be chosen to demonstrate various trajectories (see text).
% Includes Variational Equations necessary to solve Jacobian sp(4:12).
%
% Current parameters demonstrate chaotic behavior of the oscillator.
%
% IN:
% t: time (not used, but necessary for ODE solver)
% s: Input Vector
%
% OUT:
% sp: ds/dt
%
% Dimension of the system
n=3;
% Initial parameters describing chaotic behavior
mu=0.2; f=0.0; omega=1.0;
32 Chapter 1 Nonlinear Dynamics Time Series Analysis
% define equations
sp(l)=s(2);
sp(2)=mu*(l-s(l)*2)*s(2)-s(l)A3+f*cos(s(3));
sp(3)=omega;
for i=l:n
sp(n*i+l)=s(n*i+2);
sp(n*i+2)=(-2*mu*s(l)*s(2)-3*s(l)*2)*s(n*i+l)+...
mu*(l-s(l)*2)*s(n*i+2)+f*s(n*i+3);
sp(n*i+3)=0;
end
% transpose into column vector for ODE solver
sp=sp';
function [ReNormMatrix]=GramSchmidt(Matrix)
%%
Compute a set of normal-orthogonal vectors using the Gram-Schmidt algorithm.
% Matrix contains vector Matrix(:,1),Matrix(:,2),..., Matrix(:,3)
%%
IN:
% Matrix: Matrix containing original vectors.
%
% OUT:
% ReNormMatrix: Matrix with the renormalized set of vectors
Dim=size(Matrix)*[l;0];
ReNormMatrix=zeros(size(Matrix));
ReNormMatrix(:,l)=Matrix(:,1)/norm(Matrix(:,1),2);
for j=2:Dim
z=Matrix(:,j);
for i=j-l:-l:l
z=z-dot(Matrix(:,j),ReNormMatrix(:,i))*ReNormMatrix(:,i);
end
ReNormMatrix(:,j)=z/norm(z,2);
end
? Grassberger-Procaccia Algorithm
% File: CorrDim.m
%
% Implements the Grassberger-Procaccia algorithm to measure the correlation
% dimension (data assumed tobe stored in matrix 'x' f rom VanDerPolSolve . m) .
%
% Uses: Distance.m, BinFilling.m, Slope.m
% Transform x3 to sin(x3)
xnew=x;
xnew(:,3)=sin(x(:,3));
% Number of valid points in xnew matrix.
NoPoints=[size(xnew)*[l;0] size(xnew)*[l;0] size(xnew)*[l;0]];
% Calculates 32 equi-spaced bins on a logarithmic scale
[Radius]=Distance(xnew, NoPoints);
% Fills bins with the number of pairs of points with separation given by Radius .
[BinCount]=BinFilling(xnew,NoPoints,Radius);
MaxEmDim=3;
% Normalizes the matrix Radius, by its maximum value
% (for each dimension independently).
for i=l:MaxEmDim
max=Radius(32,i);
RadiusNormal(:,i)=Radius(:,i)/max;
end
% Plots the BinCount for specific Radius for all Embedding Dimensions,
figure(1);
for i-l:MaxEmDim
ifi==l
hold off
end
loglog(RadiusNormal(:,i),BinCount(:,i),'+-');
Appendix 33
if i==l
hold on
end
end
% Plot correlation Integral
ok=l;
title(strcat('Correlation Integral'));
xlabeKstrcat ( '% of Maximum Radius, MaxEmDim=',num2str (MaxEmDim) ,')'));
ylabel('% of Number of Points');
% Find the slope for all the dimensions,
for i=l:MaxEmDim
[Slopes(i), SlopeErrors(i)]=Slope(Radius(:fi),
BinCount(:,i), .6 .125);
end
% Plot Fractal Dimensions
figure(2)
Dim=l:MaxEmDim;
errorbar(Dim,Slopes,SlopeErrors, 'b*-' );
axis([O MaxEmDim+1 0 MaxEmDim]);
grid on
zoom on;
title(strcat('Data Set: ','Duffing VanDerPol Oscillator'));
xlabeK 'Dimension' ) ;
ylabel('Fractal Dimension');
function [Radius]^Distance(Portrait, NoPoints);
%
% IN:
% Portrait: Is the matrix {Number of data points, MaxEmDim} in which
% the trajectories are contained.
% NoPoints: Is the vector {1,MaxEmDim} in which the number of valid
% points for each dimension is contained.
%
% OUT:
% Radius: Is the matrix {32,MaxEmDim} in which the difference between
% the maximum and minimum distances from one point to any other , is divided
% into 32 logarithmically equal intervals (for all dimensions).
%
% Uses: DistVectPoint.m, ElimZero.m
MaxEmDim=size(Portrait)*[0;l];
NoPointsEnd=[NoPoints 1];
MinDist=ones(1,MaxEmDim)*le20;
MaxDist=zeros(1,MaxEmDim);
Radius=zeros(32,MaxEmDim);
for EmDim=l:MaxEmDim
minval=zeros(l,EmDim);
minloc=zeros(l,EmDim),
maxval=zeros(l,EmDim),
maxloc=zeros(l,EmDim),
for i=NoPointsEnd(EmDim):-1:NoPointsEnd(EmDim+1)+1
% Calculates the distances from point Portrait(i,1: EmDim) to all the
% points in matrix Portrait (1:i-1,l:EmDim)
Distances=DistVectPoint(Portrait(1:1-1,1:EmDim), Portrait(i,1:EmDim));
% Eliminates all points with distance less than Tolerance=le-10
DistanceNoZero=ElimZero(Distances, le-lO);
[minval,minloc]=min(DistanceNoZero,[],1);
[maxval,maxloc]=max(Distances,[],1);
for j=l:EmDim;
if MinDist(j)>minval(j)
MinDist(j)=minval(j);
end
if MaxDist(j)<maxval(j)
MaxDist(j)=maxval(j);
end
end
34 Chapter 1 Nonlinear Dynamics Time Series Analysis
end
end
% Fill 32 bins (equally spaced logarithmically)
for k-l:32
Radius(k,:)exp(log(MinDist)+k*(log(MaxDist)-log(MinDist))/32);
end
function [Distance]=DistVectPoint(data,point);
%
% This function calculates the distance from all elements of the matrix
% {n,MaxDim} to the point {l,MaxDim}, for all the dimensions.
%
% IN:
% data: Is a matrix {n,MaxDim} of n points to which 'point' is going
to be compared.
% point: Is a vector {1,MaxDim} that represents the 'point' in MaxDim
% dimensions.
% OUT:
% Distance: Is a matrix {n,MaxDim} that contains the distances from
% 'point' to all other points in 'data' (for dimensions 1 to MaxDim).
%
% Example: data=[ 0
% 0
% 0
% 1
% 1
% point=[ 0
% Distance=[ 0
% 0
% 0
% 1.0000
% 1.0000
Diffe=zeros(size(data));
for i=l:size(data)*[0;l]
Diffe(:,i)=data(:,i)-point(i);
end
% Calculate Euclidean distance
Diffe=Diffe.A2;
Distance=cumsum(Diffe,2);
Distance=sqrt(Distance);
0
0
1
0
1
000
1.0000
1.0000
1.4142
0
1
1
1
1
00
111
1
];
];
.0000
.4142
.4142
.7321 ]
function DistanceNoZero=ElimZero(Distance, Tolerance);
%
% Replaces all points with distance less than Tolerance with le20.
% This is necessary in the bin-filling algorithm.
%
% IN:
% Distance: Is a matrix {n,MaxDim} that contains the distances from
% 'point' to all other points in 'data' (for dimensions 1 to MaxDim).
% Tolerance: Is a scalar that determines the minimum distance to be
considered.
%
% OUT:
% DistanceNoZero: Is a matrix {n,MaxDim} equal to Distance, but with all
% elements smaller than Tolerance replaced with le20
SigDist=Distance-Tolerance;
SigDist=((sign(sign(SigDist.*-l)-0.5))+l)*le20;
DistanceNoZero=Distance+SigDist;
function [BinCount]=BinFilling(Portrait,NoPoints,Radius)
%
% IN:
% Portrait: Is the matrix {Number of data points, MaxEmDim} in which the
% trajectories are contained.
Appendix 35
% NoPoints : Is the vector {l,MaxEmDim} in which the number of points for each
% dimension is contained.
% Radius: Is the matrix {32,MaxEmDim} in which the difference between the
% maximum and minimum distances from one point to any other , is divided
into
% 32 logarithmically equal intervals (for all dimensions)
%
% OUT:
% BinCount: Is a matrix {32 ,MaxEmDim} with the total count of pair of points
% with a distance smaller than that specified by Radius for the 32
intervals.
%
% Uses: DistVectPoint.m, CountPoints.m
MaxEmDim=size(Portrait)*[0;1];
BinCount=zeros(32,MaxEmDim);
NoPointsEnd=[NoPoints 1] ;
for EmDim=l:MaxEmDim
for i=NoPointsEnd(EmDim):-1:NoPointsEnd(EmDim+1)+1
Distances=zeros(i-l,EmDim);
Distances=DistVectPoint(Portrait(1:i-1,l:EmDim),
Portrait(i,l:EmDim));
for j=l:32
BinCount(j,1:EmDim)^BinCount(j,l:EmDim)+...
CountPoints(Distances,Radius(j,l:EmDim));
end
end
end
BinCount-BinCount./(((ones(32,1)*NoPoints).*(ones(32,1)*NoPoints-l))/2);
function [CountVect]=CountPoints(Distances,Threshold);
%
% IN:
% Distance: Is amatrix {n,MaxDim} that contains the distances from 'point'
% to all other points in 'data' for dimensions 1 to MaxDim.
% Threshold: Is the upper bound on Distance.
%
% OUT:
% CountVect: Is a vector {l,MaxDim] with the count of distances smaller
% than Threshold
VectLength=length(Threshold);
NumOfPoints=size(Distances)*[l;0];
CountVect=zeros(l,VectLength);
ThresholdMatr=ones(NumOfPoints,1)*Threshold;
CountVect=sum((Distances<ThresholdMatr),1);
function [Slope, SlopeError]=Slope(RadiusV, BinCountV, center, high)
%
% This function gives the slope and error for a line (in a logarithmic scale)
% given by the points RadiusV, BinCountV. The only relevant points are the ones
% that are contained in the band center-high/2, center+high/2.
%
% The values for center define the position of the center of the band and can
% range from 0 to 1 with 1 at the top.
%
% IN:
% RadiusV: Vector with the radii limits of a specific Dimension.
% BinCountV: Vector containing the counts of pairs of elements with distance
% smaller than radius.
% Center: Center position where the slope is to be evaluated.
% High: Band size around center.
%
% OUT:
% Slope: Slope evaluated in the region center-high/2, center+high/2.
% SlopeError: Error of the evaluated fitted line to the original data.
lnRadiusV=log(RadiusV);
lnBinCountV=log(BinCountV);
36 Chapter 1 Nonlinear Dynamics Time Series Analysis
Max=O;
Min=lnBinCountV(l);
IntervalHigh=(Max-Min)*high;
Top=-((Max-Min)*(l-center))+(IntervalHigh/2);
Base=-((Max-Min)*(l-center))-(IntervalHigh/2);
k=l;
for i=l:32
if ((lnBinCountV(i)>=Base & lnBinCountV(i)<=Top))
RelDataX(k)=lnRadiusV(i);
RelDataY(k)=lnBinCountV(i);
k=k+l;
end
end
[P,S]=polyfit(RelDataX,RelDataY,l);
Slope=P(l);
SlopeError=S.normr;