sde_solver.m

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

function varargout = sde_solver(varargin)
% SDE_SOLVER provides a graphical user interface for SDE solving
% GUI-mode using:
% data=sde_solver;
% batch (no-GUI) mode using:
% data=sde_solver(f,g,agnv,awnv,init,time,solver)
% Parameters:
% f - cell array of right-side deterministic part
% g - cell array of right-side noise part
% agnv, awnv - cell arrays with additive Gaussian and additive white noise variances.
% init - - cell array of initial conditions
% time - a vector specifying the interval of integration. integration step
% must be specifying in xplicit form
% solver - solver type. now avialable only
%      1 - Milstein (Ito SDE); 
%
% data - solution matrix, first column - solution time, other columns - solution array
%
% Examples:
% Stochastic Duffing: D=sde_solver({'0.03*t*x1-x1^3'},{'0.5'},{'0'},{'0'},{'0'},'0:0.01:100','1'); 
% data=sde_solver({'-x2-x3' 'x1+0.2*x2' '0.4+x3*(x1-8)'},{'0.1' '0.1' '0.1'},{'0' '0' '0'},{'0' '0' '0'},{'1' '1' '1'},'0:0.01:100','1');
% data=sde_solver({'-x2-x3' 'x1+0.25*x2+x4' '3+x3*x1' '-0.5*x3+0.05*x4'},{'0.01*exp(-((x1+60)/10)^2-(x2/10)^2)' '0.01*exp(-((x1+60)/10)^2-(x2/10)^2)' '0.01*exp(-((x1+60)/10)^2-(x2/10)^2)' '0.01*exp(-((x1+60)/10)^2-(x2/10)^2)'},{'0' '0' '0' '0'},{'0' '0' '0' '0'},{'-15' '11' '0.2' '23'},'0:0.001:300','1');
% data=sde_solver({'10*(x2-x1)' '28*x1-x2-x1*x3' '-8/3*x3+x1*x2'},{'0.1' '0.1' '0.1'},{'0' '0' '0'},{'0' '0' '0'},{'1' '1' '1'},'0:0.01:100','1');
% See also ODE
% author: Max Logunov, [e-mail], comments welcome
% This program is a part of Lab432 software for nonlinear analysis of time
% series (freely available by request)

% Last Modified by GUIDE v2.5 14-Sep-2005 16:07:25
% last modified 8.10.06

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @sde_solver_OpeningFcn, ...
    'gui_OutputFcn',  @sde_solver_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 sde_solver_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)];
    end
    set(handles.x,'string',x);
    set(handles.f,'string',varargin{1});
    set(handles.g,'string',varargin{2});
    set(handles.agnv,'string',varargin{3});
    set(handles.annv,'string',varargin{4});
    set(handles.init,'string',varargin{5});
    set(handles.edit_length,'string',varargin{6});
    set(handles.solver,'value',varargin{7});
catch
end

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

% --- Outputs from this function are returned to the command line.
function varargout = sde_solver_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 sdeDATA

if nargout
    if handles.auto
        ok_Callback(handles.ok, [], handles);
        varargout{1} =sdeDATA;
    else
        uiwait;
        varargout{1} =sdeDATA;
    end
else
    warning('Require output argument then call sde_solver');
end

try
    close(handles.figure);
catch
end


% --- Executes on button press in ok.
function ok_Callback(hObject, eventdata, handles)
% hObject    handle to ok (see GCBO)
% handles    structure with handles and user data (see GUIDATA)
global sdeDATA
clear tmp_f
clear tmp_g
clear tmp_dg

sdeDATA=[];

x=get(handles.x,'string');
g=get(handles.g,'string'); 
init=get(handles.init,'string');
f=get(handles.f,'string');
if isempty(x)
    msgbox('No data for system');
    return
end
s = which('sde_solver');
[pathstr,name,ext,versn] = fileparts(s);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
init_val=[];
for i=1:length(x)
    init_val=[init_val; str2num(init{i})];
end

fid = fopen([pathstr '/tmp_f.m'],'w');
fid2 = fopen([pathstr '/tmp_g.m'],'w');
fprintf(fid,'function Y=tmp_f(t,X); \n');
fprintf(fid2,'function Y=tmp_g(t,X); \n');
fprintf(fid,'Y=zeros(%s,1); \n',num2str(length(x)));
fprintf(fid2,'Y=zeros(%s,1); \n',num2str(length(x)));

for i=1:length(x)
    fprintf(fid,'%s=X(%i); \n',['x' num2str(i)],i);
    fprintf(fid2,'%s=X(%i); \n',['x' num2str(i)],i);
end
for i=1:length(x)
    fprintf(fid,'Y(%s)=%s; \n',num2str(i),f{i});
    fprintf(fid2,'Y(%s)=%s; \n',num2str(i),g{i});
end
fclose(fid);
clear fid
fclose(fid2);
clear fid2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch get(handles.solver,'value')
    case 1 %Milstein method
        try
            syms dfgdsgfdf
        catch
            error('Symbolic Toolbox required. Cannot continue.');
        end
        fid = fopen([pathstr '/tmp_dg.m'],'w');
        fprintf(fid,'function Y=tmp_dg(t,X); \n');
        fprintf(fid,'Y=zeros(%s,1); \n',num2str(length(x)));
        for i=1:length(x)
            fprintf(fid,'%s=X(%i); \n',['x' num2str(i)],i);
        end
        isDgzeros=1;
        for i=1:length(x)
            r_hand=char(diff(sym(g{i}),['x' num2str(i)]));
            fprintf(fid,'Y(%s)=%s; \n',num2str(i),r_hand);
            if ~strcmp(r_hand,'0');
                isDgzeros=0;
            end
        end
        fclose(fid);
        clear fid
        time=eval(get(handles.edit_length,'string')); % time steps must be equidistant!
        dt=time(2)-time(1);
        dW=zeros(length(x),length(time));
        dW(:,2:end)=dt.^0.5.*randn(length(x),length(time)-1);
        X=inf(length(x),length(time));
        X(:,1)=init_val;
        if isDgzeros
            Vect=tmp_g(time(2),X(:,1));
            for i=2:length(time)
                X(:,i)=X(:,i-1)+dt.*tmp_f(time(i),X(:,i-1))+Vect.*dW(:,i);
                if ~isfinite(X(1,i))
                    break
                end
            end
        else            
            for i=2:length(time)
                X(:,i)=X(:,i-1)+dt.*tmp_f(time(i),X(:,i-1))+tmp_g(time(i),X(:,i-1)).*dW(:,i)+0.5.*tmp_g(time(i),X(:,i-1)).*tmp_dg(time(i),X(:,i-1)).*(dW(:,i).^2-dt);
                if ~isfinite(X(1,i))
                    break
                end
            end
        end
        clear tmp_dg
        
end
clear tmp_f
clear tmp_g


agnv=get(handles.agnv,'string');
annv=get(handles.annv,'string');
g_n=randn(size(X));
n_n=rand(size(X));
for i=1:length(x)
    g_n(i,:)=g_n(i,:)*str2num(agnv{i})^0.5;
    n_n(i,:)=(n_n(i,:)-mean(n_n(i,:)))/(var(n_n(i,:)))^.5*str2num(annv{i})^0.5;
end
X=X+g_n+n_n;

sdeDATA=[time' X'];

uiresume;


% --- 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 SDE');
if file~=0
    try
        S=load([path file]);
        set([handles.x handles.f handles.init handles.g handles.annv handles.agnv],'value',1);
        set(handles.x,'string',S.sde.dx);
        set(handles.f,'string',S.sde.f);
        set(handles.g,'string',S.sde.g);
        set(handles.agnv,'string',S.sde.agnv);
        set(handles.annv,'string',S.sde.annv);
        set(handles.init,'string',S.sde.init);
        set(handles.edit_info,'string',S.sde.info);
        set(handles.edit_length,'string',S.sde.time);
        str=get(handles.solver,'string');
        for i=1:length(str)
            if strcmp(str{i},S.sde.solver)
                set(handles.solver,'value',i);
                break
            end