lyap_exp_of_cont_sys.m

主成分分析和偏最小二乘SquaresPrincipal成分分析（ PCA ）和偏最小二乘（ PLS ）

function varargout = lyap_exp_of_cont_sys(varargin)
% LYAP_EXP_OF_CONT_SYS computes Lyapunov exponents spectra of ODE
% gui-mode using:
% [LyapExps localLyapExps]=lyap_exp_of_cont_sys;
% batch-mode using:
% [LyapExps localLyapExps]=lyap_exp_of_cont_sys(Fx,init,time,solver,tr_time,show_local)
% LyapExps=lyap_exp_of_cont_sys(Fx,init,time,solver,tr_time,show_local)
% Parameters:
% Fx - cell array of right-side ODE equations
% init - cell array of initial conditions
% time - a vector specifying the interval of integration. Must be written in form: '[t_0:step:t_end]' ('[t_0 t_end]' not allowed)
% solver - solver type:
%             1 - ode45;
%             2 - ode23;
%             3 - ode113;
%             4 - ode15s;
%             5 - ode23s;
%             6 - ode23t;
%             7 - ode23tb;
% tr_time - a vector specifying the interval of preliminary integration
% show_local -  1 - show local Lyapunov exponent on phase portrait;
%               0 - do not show local Lyapunov exponent on phase portrait; 
%
% LyapExps - Lyapunov exponents
% localLyapExps - structure with fields
%       .lle - local Lyapunov exponents
%       .time - solution time
%       .data - ODE solution array
%
% Examples:
% LyapExps=lyap_exp_of_cont_sys({'0.032*x1+0.5*x2-x3' 'x1-0.1*x2-x1^3' '0.1*x1'},{'0.01'; '-0.01' ;'-0.02'},'[0:0.02:50]',1,'[-10 0]',1);
% [LE loc]=lyap_exp_of_cont_sys({'-x2-x3' 'x1+0.2*x2' '0.4+x3*(x1-8)'},{'10'; '-4' ;'-2'},'[0:0.02:150]',1,'[-20 0]',0);
% [LE loc]=lyap_exp_of_cont_sys({'10*(x2-x1)' '28*x1-x2-x1*x3' '-8/3*x3+x1*x2'},{'1'; '2' ;'-2'},'[0:0.02:150]',1,'[-20 0]',0);

% Last Modified by GUIDE v2.5 04-Feb-2005 11:06:06
% last modified 23.03.06

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lyap_exp_of_cont_sys_OpeningFcn, ...
                   'gui_OutputFcn',  @lyap_exp_of_cont_sys_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before discrete_map_generator is made visible.
function lyap_exp_of_cont_sys_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to discrete_map_generator (see VARARGIN)

% Choose default command line output for discrete_map_generator
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);
handles.auto=0;
try
	N=length(varargin{1});
	x={};
	for i=1:N
        x{i}=['dx' num2str(i) '/dt'];
	end
	set(handles.x,'string',x);
    set(handles.f,'string',varargin{1});
    set(handles.init,'string',varargin{2});
	set(handles.edit_length,'string',varargin{3});
    set(handles.solver,'value',varargin{4});
    set(handles.e_transient,'string',varargin{5});
    set(handles.isShowMLLE,'value',varargin{6});
catch
end

if length(varargin)==6 % batch mode OK 
    handles.auto=1;
end
guidata(hObject, handles);

% --- Outputs from this function are returned to the command line.
function varargout = lyap_exp_of_cont_sys_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
global lecs locallecs

handles.isComputeLLE=nargout;
guidata(hObject, handles);

if handles.auto
    ok_Callback(handles.ok, [], handles);
    if nargout
        varargout{1} =lecs;
    end
    if nargout==2
        varargout{2} =locallecs;
    end
    close(handles.figure);
else
    if nargout
        uiwait;
        if nargout
            varargout{1} =lecs;
        end
        if nargout==2
            varargout{2} =locallecs;
        end
    end
end


% --- Executes on button press in ok.
function ok_Callback(hObject, eventdata, handles,tmp)
% hObject    handle to ok (see GCBO)
% handles    structure with handles and user data (see GUIDATA)
global lecs locallecs
lecs=[]; locallecs=[];

try 
    sym('asddaewrrjfsd');
catch
    warndlg('Symbolic toolbox not installed. Computation aborted.','Symbolic toolbox needed');
    return
end
clear tmp_sys
x=get(handles.x,'string');
init=get(handles.init,'string');
f=get(handles.f,'string');
Dim=length(f);

if isempty(x)
    msgbox('No data for system');
    return
end
s = which('lyap_exp_of_cont_sys');
[pathstr,name,ext,versn] = fileparts(s);

init_str='';
for i=1:Dim
    init_str=[init_str init{i} ' '];
end
init_str=['[' init_str(1:end-1) ']'];

fid = fopen([pathstr '/tmp_sys.m'],'w');
fprintf(fid,'function Y=tmp_sys(t,X); \n');
fprintf(fid,'Y=zeros(%s,1); \n',num2str(length(x)));
for i=1:Dim
	fprintf(fid,'%s=X(%i); \n',['x' num2str(i)],i);
end
for i=1:Dim
	fprintf(fid,'Y(%s)=%s; \n',num2str(i),f{i});
end
fclose(fid);
clear fid
 
solvers=get(handles.solver,'string');
if ~isempty(get(handles.e_transient,'string'))
	[T,Y]=eval([solvers{get(handles.solver,'value')} '(@tmp_sys,' get(handles.e_transient,'string') ',' init_str ')']); 
    Init=Y(end,:);
else
    Init=zeros(1,Dim);
	for i=1:Dim
        Init(i)=str2double(init{i});
	end
end
clear tmp_sys;

if get(handles.isShowMLLE,'value')
	InitCond='[';
	for k=1:length(Init)
        InitCond=[InitCond ' ' num2str(Init(k))];
	end
	InitCond=[InitCond ']'];
	[Tplot,Yplot]=eval([solvers{get(handles.solver,'value')} '(@tmp_sys,' get(handles.edit_length,'string') ',' InitCond ')']); 
end
clear tmp_sys;

F=[]; X=[];
for i=1:Dim
    F=[F; sym(f{i})];
    X=[X sym(['x' num2str(i)])];
end
J=jacobian(F,X);

Xer=[];
for i=1:Dim^2
    Xer=[Xer; sym(['x' num2str(i) 'er'])];
end

system=[];
for i=1:Dim
    system=[system; J*Xer(1+(i-1)*Dim:i*Dim)];
end
%delete tmp_sys.m
fid = fopen([pathstr '/tmp_sys.m'],'a');
fprintf(fid,'Y=[Y; zeros(%s,1)]; \n',num2str(length(system)));
for i=1:length(Xer)
    fprintf(fid,'%s=X(%i); \n',char(Xer(i)),Dim+i);
end
for i=1:length(system)
    fprintf(fid,'Y(%s)=%s; \n',num2str(Dim+i),char(system(i)));
end

fclose(fid);
clear fid
clear tmp_sys

Vect=eye(Dim);
S=zeros(Dim,1);
Y=Init';
Y=[Y; Vect(1:end)']';

M=str2num(get(handles.edit_length,'string')); %#ok<ST2NM>
dt=M(2)-M(1);

isCompLLE=handles.isComputeLLE|get(handles.isShowMLLE,'value');
if isCompLLE
	locallecs=zeros(length(1:2:length(M)-2),Dim);
    time=zeros(length(1:2:length(M)-2),1);
    sol=zeros(length(1:2:length(M)-2),Dim);
end
counter=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:2:length(M)-2
%     wwwhan=waitbar(i/length(M));
	InitCond=['[' num2str(Y(:)','%20.15f') ']'];
    
    [T,Y]=eval([solvers{get(handles.solver,'value')} '(@tmp_sys,[' num2str(M(i),8) ',' num2str(M(i+1),8) ',' num2str(M(i+2),8) '],' InitCond ',1e-5)']); 
	if isCompLLE
        sol(counter,:)=Y(2,1:Dim); %%ATTENTION!! -this line is correct only if length(T)==3!!
	end
    Y=Y(end,:);
    vectors=Y(end-Dim^2+1:end);
    for t=1:Dim
        Vect(:,t)=vectors(1+(t-1)*Dim:t*Dim)';
    end
    localExp=log(sum(Vect(:,1).^2)^.5);
    S(1)=S(1)+localExp;
	if isCompLLE
        locallecs(counter,1)=localExp/dt; 
        time(counter)=T(1)+(T(end)-T(1))/2;
	end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     if log(sum(Vect(:,1).^2)^.5)>0
%         plot3(Y(1),Y(2),Y(3),'b')
%         grid on
%         hold on 
%         plot3([Y(1), Y(1)+Vect(1,1)],[Y(2), Y(2)+Vect(2,1)],[Y(3), Y(3)+Vect(3,1)],'r')
%         drawnow
%     end
%     colors={'y' 'g' 'b' 'r'};
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Vect(:,1)=Vect(:,1)./sum(Vect(:,1).^2)^.5;
    for j=2:Dim
        tmp_v=zeros(Dim,1);
        for k=1:j-1
            tmp_v=tmp_v+sum(Vect(:,j).*Vect(:,k)).*Vect(:,k);
        end
        Vect(:,j)=Vect(:,j)-tmp_v;
        localExp=log(sum(Vect(:,j).^2)^.5);
        S(j)=S(j)+localExp;
		if isCompLLE
            locallecs(counter,j)=localExp/dt;
		end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     if log(sum(Vect(:,j).^2)^.5)>0
%         plot3([Y(1), Y(1)+Vect(1,j)],[Y(2), Y(2)+Vect(2,j)],[Y(3), Y(3)+Vect(3,j)],colors{j})
%         drawnow
%     end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Vect(:,j)=Vect(:,j)./sum(Vect(:,j).^2)^.5;
    end
    counter=counter+1;
    Y=[Y(1:Dim) Vect(1:end)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isCompLLE
	locallecs.lle=locallecs/2;
	locallecs.time=time;
    locallecs.sol=sol;
end
% delete(wwwhan);
clear tmp_sys
fprintf('\nLyapunov exponents:\n ');
Le=S./length(M)/dt;
for i=1:Dim
    fprintf('       %f\n',Le(i));
end
lecs=Le;
if get(handles.isShowMLLE,'value')
    F=figure('Name','Maximal local Lyapunov exponent' ,'NumberTitle','off','color',[1 1 1]);
    Tplot=Tplot(2:2:end);
    Yplot=Yplot(2:2:end,:);
    L=min([length(Tplot) length(locallecs.lle(:,1))]);
    Tplot=Tplot(1:L); Yplot=Yplot(1:L,:); locallecs.lle=locallecs.lle(1:L,:);
    [i j]=size(Yplot);
    if j==1
        h = patch([Tplot' fliplr(Tplot')], [Yplot' fliplr(Yplot')], [locallecs.lle' fliplr(locallecs.lle')]); grid on;