region_of_stability

 function to analyse the region of stability in user defined boundaries

 input:   o SETTINGS....struct with all informations about the current
                        variables

 output:  o none

 notes:   o variables from the last run are saved as file for next run

 authors: o Christian Jäkel (University of Technology Dresden)
          o Daniel Klawitter (University of Technology Dresden)

0001 function region_of_stability(SETTINGS)
0002 % function to analyse the region of stability in user defined boundaries
0003 %
0004 % input:   o SETTINGS....struct with all informations about the current
0005 %                        variables
0006 %
0007 % output:  o none
0008 %
0009 % notes:   o variables from the last run are saved as file for next run
0010 %
0011 % authors: o Christian Jäkel (University of Technology Dresden)
0012 %          o Daniel Klawitter (University of Technology Dresden)
0013 
0014 promt={'lower intervall boundary on the x-axis           ','upper intervall boundary on the x-axis',...
0015     'lower intervall boundary on the y-axis','upper intervall boundary on the y-axis',...
0016     'fineness for the x-axis','fineness for the y-axis','size of the dots  (higher number --> bigger dots)'};
0017 if exist('data\default_value_stability.mat','file')
0018     
0019     load data\default_value_stability.mat
0020 else
0021     a1=-5;
0022     a2=2;
0023     a3=-4;
0024     a4=4;
0025     a5=0.05;
0026     a6=0.05;
0027     a7=0.5;
0028     
0029 end
0030 
0031 defAns_cell={num2str(a1),num2str(a2),...
0032     num2str(a3),num2str(a4)...
0033     num2str(a5),num2str(a6),num2str(a7)};
0034 
0035 get_input=inputdlg(promt,'parameters',1,defAns_cell);
0036 
0037 a1=str2num(get_input{1});
0038 a2=str2num(get_input{2});
0039 a3=str2num(get_input{3});
0040 a4=str2num(get_input{4});
0041 a5=str2num(get_input{5});
0042 a6=str2num(get_input{6});
0043 a7=str2num(get_input{7});
0044 
0045 % save default values for next run
0046 save data\default_value_stability.mat a1 a2 a3 a4 a5 a6 a7
0047 
0048 if SETTINGS.isRKV
0049     B =SETTINGS.B;
0050     c =SETTINGS.c;
0051     a =SETTINGS.a;
0052     s = size(B,1);
0053 else
0054     rho    = SETTINGS.rho;
0055     sigma  = SETTINGS.sigma;
0056 end
0057 
0058 % the region witch is tested for stability
0059 X = a1:a5:a2;
0060 Y = a3:a6:a4;
0061 i = sqrt(-1);
0062 
0063 count_idx = 0;
0064 if SETTINGS.isRKV
0065     % Runge-Kutta-Method
0066     for x_idx = 1:length(X)
0067         for y_idx = 1:length(Y)
0068             if abs(1+c*(eye(s)-(X(x_idx)+i*Y(y_idx))*B)^(-1)*ones(1,s)'*(X(x_idx)+i*Y(y_idx)))<=1  % the stability function of a runge-kutta-method in general form is
0069                 count_idx             = count_idx+1;                                               % eveluatet at the specivic grid points.
0070                 plot_vec_x(count_idx) = X(x_idx);
0071                 plot_vec_y(count_idx) = Y(y_idx);
0072             end
0073         end
0074     end
0075 else
0076     % linear multistep method
0077     for x_idx = 1:length(X)
0078         for y_idx = 1:length(Y)
0079             Wurzeln = roots(rho-(X(x_idx)+i*Y(y_idx))*sigma);
0080             if abs(Wurzeln) <= 1 & (length(find(Wurzeln==1))==1 | length(find(Wurzeln==1))==0)% root condition
0081                 count_idx             = count_idx+1;
0082                 plot_vec_x(count_idx) = X(x_idx);
0083                 plot_vec_y(count_idx) = Y(y_idx);
0084             end
0085         end
0086     end
0087 end
0088 
0089 if ~exist('plot_vec_y','var')
0090     errordlg('For these boundaries the region of stability is empty!','error')
0091 else
0092     figure
0093     plot(plot_vec_x,plot_vec_y,'*','MarkerSize',a7)
0094 end
0095 end