www.gusucode.com > Matlab动力系统和时间序列分析工具箱 > Matlab动力系统和时间序列分析工具箱/lab432/common/lyap_exp_of_cont_sys.m
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;
set(h, 'EdgeColor', 'Flat');
xlabel('t');
ylabel('x1');
end
if j==2
h = patch([Yplot(:,1)' fliplr(Yplot(:,1)')], [Yplot(:,2)' fliplr(Yplot(:,2)')], [max(locallecs.lle') fliplr(max(locallecs.lle'))]); grid on;
set(h, 'EdgeColor', 'Flat');
xlabel('x1');
ylabel('x2');
end
if j>2
h = patch([Yplot(:,1)' fliplr(Yplot(:,1)')], [Yplot(:,2)' fliplr(Yplot(:,2)')], [Yplot(:,3)' fliplr(Yplot(:,3)')], [max(locallecs.lle') fliplr(max(locallecs.lle'))]); grid on;
set(h, 'EdgeColor', 'Flat');
xlabel('x1');
ylabel('x2');
zlabel('x3');
view(45,45);
end
end
% --- Executes on button press in cancel.
function cancel_Callback(hObject, eventdata, handles)
% hObject handle to cancel (see GCBO)
% handles structure with handles and user data (see GUIDATA)
close(handles.figure);
% --- Executes during object creation, after setting all properties.
function edit_length_CreateFcn(hObject, eventdata, handles)
% hObject handle to edit_length (see GCBO)
% handles empty - handles not created until after all CreateFcns called
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
function edit_length_Callback(hObject, eventdata, handles)
% hObject handle to edit_length (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% --- Executes on button press in load.
function load_Callback(hObject, eventdata, handles)
% hObject handle to load (see GCBO)
% handles structure with handles and user data (see GUIDATA)
[file,path] = uigetfile2('*.mat','Load System');
if file~=0
try
S=load([path file]);
set([handles.x handles.f handles.init],'value',1);
set(handles.x,'string',S.sys.dx);
set(handles.f,'string',S.sys.f);
set(handles.init,'string',S.sys.init);
set(handles.edit_info,'string',S.sys.info);
set(handles.edit_length,'string',S.sys.time);
set(handles.e_transient,'string',S.sys.tr_time);
str=get(handles.solver,'string');
for i=1:length(str)
if strcmp(str{i},S.sys.solver)
set(handles.solver,'value',i);
break
end
end
catch
errordlg('Wrong file format','File error');
end
end
% --- Executes on button press in save.
function save_Callback(hObject, eventdata, handles)
% hObject handle to save (see GCBO)
% handles structure with handles and user data (see GUIDATA)
[file,path] = uiputfile('*.mat','Save System As');
if file~=0
sys.dx=get(handles.x,'string');
sys.f=get(handles.f,'string');
tmp=sys.dx;
for i=1:length(tmp)
tmp{i}='0';
end
sys.ignv=tmp;
sys.innv=tmp;
sys.agnv=tmp;
sys.annv=tmp;
sys.init=get(handles.init,'string');
sys.info=get(handles.edit_info,'string');
sys.tr_time=get(handles.e_transient,'string');
sys.time=get(handles.edit_length,'string');
str=get(handles.solver,'string');
sys.solver=str{get(handles.solver,'value')};
sys.addit_param='';
save([path file],'sys');
end
% --- Executes during object creation, after setting all properties.
function edit_info_CreateFcn(hObject, eventdata, handles)
% hObject handle to edit_info (see GCBO)
% handles empty - handles not created until after all CreateFcns called
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
function edit_info_Callback(hObject, eventdata, handles)
% hObject handle to edit_info (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% --- Executes during object creation, after setting all properties.
function x_CreateFcn(hObject, eventdata, handles)
% hObject handle to x (see GCBO)
% handles empty - handles not created until after all CreateFcns called
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on selection change in x.
function x_Callback(hObject, eventdata, handles)
% hObject handle to x (see GCBO)
% handles structure with handles and user data (see GUIDATA)
set([handles.x handles.f handles.init],'value',get(hObject,'Value'));
% --- Executes during object creation, after setting all properties.
function f_CreateFcn(hObject, eventdata, handles)
% hObject handle to f (see GCBO)
% handles empty - handles not created until after all CreateFcns called
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes during object creation, after setting all properties.
function init_CreateFcn(hObject, eventdata, handles)
% hObject handle to init (see GCBO)
% handles empty - handles not created until after all CreateFcns called
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on button press in add.
function add_Callback(hObject, eventdata, handles)
% hObject handle to add (see GCBO)
% handles structure with handles and user data (see GUIDATA)
N=length(get(handles.x,'string'));
x={};
for i=1:N
x{i}=['dx' num2str(i) '/dt'];
end
set(handles.x,'string',x);
prompt = {['Enter equation dx' num2str(N+1) '/dt=...'],...
'Enter initial condition:'};
dlg_title = 'Add new equation';
num_lines= 1;
def = {'','0.5'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
if ~isempty(answer)
x=get(handles.x,'string');
f=get(handles.f,'string');
init=get(handles.init,'string');
x{N+1}=['dx' num2str(N+1) '/dt'];
f{N+1}=answer{1};
init{N+1}=answer{2};
set([handles.x, handles.f, handles.init],'value',N+1);
set(handles.x,'string',x);
set(handles.f,'string',f);
set(handles.init,'string',init);
end
% --- Executes on button press in modify.
function modify_Callback(hObject, eventdata, handles)
% hObject handle to modify (see GCBO)
% handles structure with handles and user data (see GUIDATA)
x=get(handles.x,'string');
f=get(handles.f,'string');
init=get(handles.init,'string');
selection=get(handles.x,'value');
if isempty(x)
x{1}='dx1/dt'; f{1}=''; init{1}='0.7';
end
prompt = {['Enter equation ' x{selection} '=...'],...
'Enter initial condition:'};
dlg_title = 'Modify equation';
num_lines= 1;
def = {f{selection},init{selection}};
answer = inputdlg(prompt,dlg_title,num_lines,def);
if ~isempty(answer)
f{selection}=answer{1};
init{selection}=answer{2};
set(handles.x,'string',x);
set(handles.f,'string',f);
set(handles.init,'string',init);
end
% --- Executes on button press in delete.
function delete_Callback(hObject, eventdata, handles)
% hObject handle to delete (see GCBO)
% handles structure with handles and user data (see GUIDATA)
x=get(handles.x,'string');
f=get(handles.f,'string');
init=get(handles.init,'string');
selection=get(handles.x,'value');
if isempty(x)
return
end
n_f={}; n_init={}; n_x={};
k=1;
for i=1:length(x)
if i~=selection
n_f{k}=f{i};
n_init{k}=init{i};
n_x{k}=x{i};
k=k+1;
end
end
set([handles.x, handles.f, handles.init],'value',max([1 selection-1]));
set(handles.x,'string',n_x);
set(handles.f,'string',n_f);
set(handles.init,'string',n_init);
% --- Executes during object creation, after setting all properties.
function solver_CreateFcn(hObject, eventdata, handles)
% hObject handle to solver (see GCBO)
% handles empty - handles not created until after all CreateFcns called
% Hint: popupmenu controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject,'string',{'ode45' 'ode23' 'ode113' 'ode15s' 'ode23s' 'ode23t' 'ode23tb'});
% --- Executes on selection change in solver.
function solver_Callback(hObject, eventdata, handles)
% hObject handle to solver (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hints: contents = get(hObject,'String') returns solver contents as cell array
% contents{get(hObject,'Value')} returns selected item from solver
% --- Executes during object creation, after setting all properties.
function e_transient_CreateFcn(hObject, eventdata, handles)
% hObject handle to e_transient (see GCBO)
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
function e_transient_Callback(hObject, eventdata, handles)
% hObject handle to e_transient (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of e_transient as text
% str2double(get(hObject,'String')) returns contents of e_transient as a double
% --- Executes on button press in isComputeLLE.
function isComputeLLE_Callback(hObject, eventdata, handles)
% hObject handle to isComputeLLE (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hint: get(hObject,'Value') returns toggle state of isComputeLLE
% --- Executes on button press in isShowMLLE.
function isShowMLLE_Callback(hObject, eventdata, handles)
% hObject handle to isShowMLLE (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hint: get(hObject,'Value') returns toggle state of isShowMLLE