odeoutp.m

关于混沌系统的李氏指数计算等混沌系统中重要参数计算的代码

                  rr = eval(DS(1).poincare_eq);
                  rr1=rr;
                  XK=X; TK=t;
                  XL=TRJ_bufer(bufer_i-1,1:neq);
                  TL=Time_bufer(bufer_i-1);
                  if (rr*DS(1).poincare_cur)<0
                      k = 0;
                      while (abs(rr)>DS(1).rel_error) & (k<=10)
                          k = k+1;
                          X = oderhs(TK,XK,P);
                          t = TK;
                          df_plane=eval([DS(1).poincare_eq '+' DS(1).poincare_map{2}]);
                          TK1 = TK - rr/df_plane; 
                          rr = TK1 - TL;
                          Poinc_options = odeset('RelTol',DS(1).rel_error,'AbsTol',DS(1).abs_error,'MaxStep',DS(1).max_step,...
                              'Refine',2,'InitialStep',rr);                          
                          [TT,XT]=integrator(DS(1).method_int,@oderhs,[TL TK1],XL,Poinc_options,P);
                          kXT = size(XT);
                          t = TT(kXT(1));
                          X = XT(kXT(1),1:neq);
                          rr = eval(DS(1).poincare_eq);
                          XK = X;
                          TK = t;
                          if (rr*DS(1).poincare_cur)>0
                              XL = X;
                              TL = t;
                          end
                      end;
% Output if all is OK                      
                      if k<10 
                          DS(1).poincare_cur = rr1; 
                          DS(1).poincare_nmbr = DS(1).poincare_nmbr + 1;
                          nwind = size(DS.windows); 
                          nwind = nwind(2)-1; 
                          for i=1:nwind 
                              ff = get(DS(1).windows(i+1),'UserData');
                              if ff.poinc_out==1
                                  if ff.type == 1                         
                                      Xp = eval(ff.Xvalue); 
                                      Yp = eval(ff.Yvalue);
                                      if DS(1).poincare_nmbr==1 
                                          set(ff.poinc_trj,'Xdata',Xp,'Ydata',Yp);
                                      else 
                                          ux=get(ff.poinc_trj,'XData');
                                          uy=get(ff.poinc_trj,'YData');
                                          set(ff.poinc_trj,'Xdata',[ux Xp],'Ydata',[uy Yp]);
                                      end 
                                      drawnow; 
                                  else
                                      Xp = eval(ff.Xvalue); 
                                      Yp = eval(ff.Yvalue);
                                      Zp = eval(ff.Zvalue);
                                      if DS(1).poincare_nmbr==1 
                                          set(ff.poinc_trj,'Xdata',Xp,'Ydata',Yp,'Zdata',Zp);
                                      else 
                                          ux=get(ff.poinc_trj,'XData');
                                          uy=get(ff.poinc_trj,'YData');
                                          uz=get(ff.poinc_trj,'ZData');
                                          set(ff.poinc_trj,'Xdata',[ux Xp],'Ydata',[uy Yp],'Zdata',[uz Zp]);
                                      end 
                                  end
                              end
                          end
                      end;
                  end
              case 3
% Poincare map by secant to positive direction                  
                  X = Xsol(1:neq,1)';
                  t = time(1);
                  rr = eval(DS(1).poincare_eq);
                  rr1=rr;
                  XK=X; TK=t;
                  XL=TRJ_bufer(bufer_i-1,1:neq);
                  TL=Time_bufer(bufer_i-1);
                  if ( (rr*DS(1).poincare_cur)<0 ) & (rr>0)
                      k = 0;
                      while (abs(rr)>DS(1).rel_error) & (k<=10)
                          k = k+1;
                          X = oderhs(TK,XK,P);
                          t = TK;
                          df_plane=eval([DS(1).poincare_eq '+' DS(1).poincare_map{2}]);
                          TK1 = TK - rr/df_plane; 
                          rr = TK1 - TL;
                          Poinc_options = odeset('RelTol',DS(1).rel_error,'AbsTol',DS(1).abs_error,'MaxStep',DS(1).max_step,...
                              'Refine',2,'InitialStep',rr);                          
                          [TT,XT]=integrator(DS(1).method_int,@oderhs,[TL TK1],XL,Poinc_options,P);
                          kXT = size(XT);
                          t = TT(kXT(1));
                          X = XT(kXT(1),1:neq);
                          rr = eval(DS(1).poincare_eq);
                          XK = X;
                          TK = t;
                          if (rr*DS(1).poincare_cur)>0
                              XL = X;
                              TL = t;
                          end
                      end;
% Output if all is OK                      
                      if k<10 
                          DS(1).poincare_cur = rr1; 
                          DS(1).poincare_nmbr = DS(1).poincare_nmbr + 1;
                          nwind = size(DS.windows); 
                          nwind = nwind(2)-1; 
                          for i=1:nwind 
                              ff = get(DS(1).windows(i+1),'UserData');
                              if ff.poinc_out==1
                                  if ff.type == 1                         
                                      Xp = eval(ff.Xvalue); 
                                      Yp = eval(ff.Yvalue);
                                      if DS(1).poincare_nmbr==1 
                                          set(ff.poinc_trj,'Xdata',Xp,'Ydata',Yp);
                                      else 
                                          ux=get(ff.poinc_trj,'XData');
                                          uy=get(ff.poinc_trj,'YData');
                                          set(ff.poinc_trj,'Xdata',[ux Xp],'Ydata',[uy Yp]);
                                      end 
                                      drawnow; 
                                  else
                                      Xp = eval(ff.Xvalue); 
                                      Yp = eval(ff.Yvalue);
                                      Zp = eval(ff.Zvalue);
                                      if DS(1).poincare_nmbr==1 
                                          set(ff.poinc_trj,'Xdata',Xp,'Ydata',Yp,'Zdata',Zp);
                                      else 
                                          ux=get(ff.poinc_trj,'XData');
                                          uy=get(ff.poinc_trj,'YData');
                                          uz=get(ff.poinc_trj,'ZData');
                                          set(ff.poinc_trj,'Xdata',[ux Xp],'Ydata',[uy Yp],'Zdata',[uz Zp]);
                                      end 
                                  end
                              end
                          end
                      end;
                  end
              case 4
% Poincare map by secant to negative direction                  
                  X = Xsol(1:neq,1)';
                  t = time(1);
                  rr = eval(DS(1).poincare_eq);
                  rr1=rr;
                  XK=X; TK=t;
                  XL=TRJ_bufer(bufer_i-1,1:neq);
                  TL=Time_bufer(bufer_i-1);
                  if ( (rr*DS(1).poincare_cur)<0 ) & ( rr<0 )
                      k = 0;
                      while (abs(rr)>DS(1).rel_error) & (k<=10)
                          k = k+1;
                          X = oderhs(TK,XK,P);
                          t = TK;
                          df_plane=eval([DS(1).poincare_eq '+' DS(1).poincare_map{2}]);
                          TK1 = TK - rr/df_plane; 
                          rr = TK1 - TL;
                          Poinc_options = odeset('RelTol',DS(1).rel_error,'AbsTol',DS(1).abs_error,'MaxStep',DS(1).max_step,...
                              'Refine',2,'InitialStep',rr);                          
                          [TT,XT]=integrator(DS(1).method_int,@oderhs,[TL TK1],XL,Poinc_options,P);
                          kXT = size(XT);
                          t = TT(kXT(1));
                          X = XT(kXT(1),1:neq);
                          rr = eval(DS(1).poincare_eq);
                          XK = X;
                          TK = t;
                          if (rr*DS(1).poincare_cur)>0
                              XL = X;
                              TL = t;
                          end
                      end;
% Output if all is OK                      
                      if k<10 
                          DS(1).poincare_cur = rr1; 
                          DS(1).poincare_nmbr = DS(1).poincare_nmbr + 1;
                          nwind = size(DS.windows); 
                          nwind = nwind(2)-1; 
                          for i=1:nwind 
                              ff = get(DS(1).windows(i+1),'UserData');
                              if ff.poinc_out==1
                                  if ff.type == 1                         
                                      Xp = eval(ff.Xvalue); 
                                      Yp = eval(ff.Yvalue);
                                      if DS(1).poincare_nmbr==1 
                                          set(ff.poinc_trj,'Xdata',Xp,'Ydata',Yp);
                                      else 
                                          ux=get(ff.poinc_trj,'XData');
                                          uy=get(ff.poinc_trj,'YData');
                                          set(ff.poinc_trj,'Xdata',[ux Xp],'Ydata',[uy Yp]);
                                      end 
                                      drawnow; 
                                  else
                                      Xp = eval(ff.Xvalue); 
                                      Yp = eval(ff.Yvalue);
                                      Zp = eval(ff.Zvalue);
                                      if DS(1).poincare_nmbr==1 
                                          set(ff.poinc_trj,'Xdata',Xp,'Ydata',Yp,'Zdata',Zp);
                                      else 
                                          ux=get(ff.poinc_trj,'XData');
                                          uy=get(ff.poinc_trj,'YData');
                                          uz=get(ff.poinc_trj,'ZData');
                                          set(ff.poinc_trj,'Xdata',[ux Xp],'Ydata',[uy Yp],'Zdata',[uz Zp]);
                                      end 
                                  end
                              end
                          end
                      end;
                  end
          end
      end
      
      X = Xsol(1:neq,1)';
      t = time(1);
      if DS(1).poincare_do>1
          DS(1).poincare_cur=eval(DS(1).poincare_eq);
      end

end;