Диссертация (1144294), страница 18
Текст из файла (страница 18)
png ’ ]) ;close ( fig ) ;2324endЛистинг В.7:experiment1.m – функция для моделирования эксперимента№1.1function experiment123r = 0; sigma = 35; eps = 0.001; d = -35 - eps ; b = 70 - eps ;4567function OUT = J ( X )OUT = [ - sigma sigma 0 ; r - X (3) -d -X (1) ; X (2) X (1) -b ];end89S0 = [0 , 0 , 0];126101112S12_XY = sqrt ( b *( r - d ) ) ;S12_Z = r - d ;131415[V , D ] = eig ( J ( S0 ) ) ;[ ~ , IX ] = sort ( diag (D ) , ’ descend ’) ;161718S1 = [ S12_XY , S12_XY , S12_Z ];S2 = [ - S12_XY , - S12_XY , S12_Z ];192021delta = 0.01;x_unst_01 = S0 + delta *V (: , IX (1) ) ’;2223epsTolArr = [1 e -8 , 1e -10 , 1e -12];2425saveDir = ’ ./ EXP1 ’;26272829if ~ exist ( saveDir , ’ dir ’)mkdir ( saveDir ) ;end3031for iCurrTol = 1 : numel ( epsTolArr )323334rel_tol = epsTolArr ( iCurrTol ) ; abs_tol = epsTolArr ( iCurrTol ) ;options = odeset ( ’ RelTol ’ , rel_tol , ’ AbsTol ’ , abs_tol ) ;3536T = 200; currT = T;3738currStartPoint = x_unst_01 ; calculetedTraj = []; calcT = [];3940414243saveDirTol = [ saveDir , ’/ tol = ’ ,%0.1 e ’) ];if ~ exist ( saveDirTol , ’ dir ’)mkdir ( saveDirTol ) ;endnum2str ( epsTolArr ( iCurrTol ) , ’44454647for iEndTime = 1 : 3[t , x_ust1 ] = ode45 (@ (t , x ) glorenzSyst (t , x , sigma , r , b, d ) , [0 , T ] , currStartPoint , options ) ;127calculetedTraj = [ calculetedTraj ; x_ust1 ];calcT = [ calcT ; t ];48495051525354555657585960616263646566h = figure ( ’ Visible ’ , ’ off ’) ; hold on ;xStart = plot3 ( x_unst_01 (1) , x_unst_01 (2) , x_unst_01 (3) ,’.
’ , ’ markersize ’ , 15 , ’ Color ’, ’ black ’) ;hS0 = plot3 ( S0 (1) , S0 (2) , S0 (3) , ’. ’ , ’ markersize ’ , 15 , ’Color ’ , ’ red ’) ;text ( S0 (1) , S0 (2) , S0 (3) , ’ S_0 ’ , ’ fontsize ’ , 18) ;hS1 = plot3 ( S1 (1) , S1 (2) , S1 (3) , ’. ’ , ’ markersize ’ , 15 , ’Color ’ , ’ red ’) ;text ( S1 (1) , S1 (2) , S1 (3) , ’ S_1 ’ , ’ fontsize ’ , 18) ;hold on ;hS2 = plot3 ( S2 (1) , S2 (2) , S2 (3) , ’.
’ , ’ markersize ’ , 15 , ’Color ’ , ’ red ’) ;text ( S2 (1) , S2 (2) , S2 (3) , ’ S_2 ’ , ’ fontsize ’ , 18) ;c = 1: numel ( calcT ) ;% # colorshAttr = surf ([ calculetedTraj (: ,1) , calculetedTraj (: ,1) ] ,[ calculetedTraj (: ,2) , calculetedTraj (: ,2) ],[ calculetedTraj (: ,3) , calculetedTraj (: ,3) ],[ c (:) , c (:) ] , ’ EdgeColor ’ , ’ flat ’ , ’ FaceColor ’ , ’ none ’) ;% http :// life - prog . ru / view_zam . php ? id =54colormap ( jet ( numel ( calcT ) ) ) ;67686970grid on ; axis auto ;xlabel ( ’x ’) ; ylabel ( ’y ’); zlabel ( ’z ’) ;view ([30 ,15]) , title ([ ’ tol = ’ , num2str ( epsTolArr (iCurrTol ) , ’ %0.1 e ’) , ’; T = ’ , num2str ( currT ) ]) ;71saveas (h , [ saveDirTol , ’/ T = ’ , num2str ( currT ) ’. png ’ ]) ;7273% UpdatecurrStartPoint = x_ust1 ( end , :) ;currT = currT + T ;747576end7778end7980end128Листинг В.8:experiment2.m – функция для моделирования эксперимента№2.1function experiment_chaos23global sigma r b d45r = 0; sigma = 35; eps = 1e -4; d = -35 + eps ; b = 70 + eps ; % exp 16789function OUT = J ( X )OUT = [ - sigma sigma 0 ; r - X (3) -d -X (1) ; X (2) X (1) -b ];end1011S0 = [0 , 0 , 0];121314S12_XY = sqrt ( b *( r - d ) ) ;S12_Z = r - d ;151617[V , D ] = eig ( J ( S0 ) )[ ~ , IX ] = sort ( diag (D ) , ’ descend ’) ;181920S1 = [ S12_XY , S12_XY , S12_Z ];S2 = [ - S12_XY , - S12_XY , S12_Z ];212223[ V1 , D1 ] = eig ( J ( S1 ))[ ~ , IX1 ] = sort ( real ( diag ( D1 ) ) , ’ descend ’) ;2425y0 = 10; x0 = 0; z0 = 0; x_unst_01 = [ x0 , y0 , z0 ];262728eps_tol = 1e -8; rel_tol = eps_tol ; abs_tol = eps_tol ;options = odeset ( ’ RelTol ’ , rel_tol , ’ AbsTol ’ , abs_tol );293031T = 500;[t , x_ust0 ] = ode45 ( @glorenzSyst , [0 , T ], x_unst_01 , options ) ;323334x_unst_01 = x_ust0 ( end ,:) ;[t , x_ust0 ] = ode45 ( @glorenzSyst , [0 , T ], x_unst_01 , options ) ;353637f1 = figure (1) ; hold on ;hS0 = plot3 ( S0 (1) , S0 (2) , S0 (3) , ’.
’ , ’ markersize ’ , 15 , ’ Color ’ , ’red ’) ;1293839hS1 = plot3 ( S1 (1) , S1 (2) , S1 (3) , ’. ’ , ’ markersize ’ , 15 , ’ Color ’ , ’red ’) ;4041hS2 = plot3 ( S2 (1) , S2 (2) , S2 (3) , ’. ’ , ’ markersize ’ , 15 , ’ Color ’ , ’red ’) ;4243444546474849c = 1: numel ( t ) ;h = surf ([ x_ust0 (: ,1) , x_ust0 (: ,1) ] , [ x_ust0 (: ,2) , x_ust0 (: ,2) ] , [x_ust0 (: ,3) , x_ust0 (: ,3) ] , [c (:) , c (:) ] , ’ EdgeColor ’ , ’ flat ’ , ’FaceColor ’ , ’ none ’) ;colormap ( jet ( numel ( t) ) );grid on ; axis auto ;xlabel ( ’x ’) ; ylabel ( ’y ’) ; zlabel ( ’z ’) ;hold off ;view ([20 ,50]) ;5051endВ.2Реализация алгоритма моделирования и анализа поведения сепаратрис седла и поведения отображенийПуанкаре системы лоренцевского типаВ.2.1 Алгоритм численного сканирования области параметров ианализа поведения сепаратрис седловых состояний равновесияЛистинг В.9:runDeltaBetaGridScan.m – скрипт для запуска численнойпроцедуры сканирования области параметров и исследования поведениясепаратрис седловых состояний равновесия.1clear global ;23DIR = ’ ./ GRID_TEST ’;456delta = 0.5;13078betaStep = 0.01;betaSpan = betaStep : betaStep : 2 + delta - betaStep ;91011sStep = 0.01;sSpan = sStep : sStep : 1 - sStep ;12131415% Choosing the vicinity for points near the equilibria S0 , S1 , S2eqEpsVic .
S0 = 1e -16;eqEpsVic . S12 = 1e -2;1617acc = 1e -16; rel_tol = acc ; abs_tol = acc ; InitialStep = acc /10;18192021222324radiusInf = 1000; S1 = [1 ,0 ,0]; epsVicS1 = 1e -1;ODE . solverName = ’ ode45 ’;ODE . options = odeset ( ’ RelTol ’ , rel_tol , ’ AbsTol ’ , abs_tol , ...’ InitialStep ’ , InitialStep , ’ NormControl ’ , ’ on ’ , ...’ Events ’ , @ (t , x ) trajBehaviorEvents (t , x , radiusInf , S1 ,epsVicS1 ) ) ;ODE . lim = true ;252627% Integration time :tEnd .
trans = 4000; tEnd . lim = 1000;282930% Log options :logOption = [];3132333435363738% Plot options :plotOptions . symmetric = false ;plotOptions . color = true ;plotOptions . figVisible = ’ off ’;plotOptions . figSave = true ;plotOptions . axis = [];plotOptions . view = [];394041% Main routing :deltaBetaGridScan ( DIR , delta , betaSpan , sSpan , eqEpsVic , ODE , tEnd ,logOption , plotOptions ) ;Листинг В.10:trajBehaviorEvents.m – функция, для выявления событий”траектория вне шара” и ”траектория в -окрестности точки”.1311234567function [ value , isterminal , direction ] = trajBehaviorEvents ( ~ , x ,radiusInf , xEq , epsEq )value = [ sqrt ( x (1) ^2 + x (2) ^2 + x (3) ^2) - radiusInf ;sqrt (( x (1) - xEq (1) ) ^2 + ...( x (2) - xEq (2) ) ^2 + ...( x (3) - xEq (3) ) ^2) - epsEq ];isterminal = [1; 1];direction = [0; 0];Листинг В.11:deltaBetaGridScan.m – численная процедура длясканирования области параметров (, ) и исследования поведения сепаратрисседловых состояний равновесия.1function deltaBetaGridScan ( DIR , delta , betaList , sList , eqEpsVic ,ODE , tEnd , logOption , plotOptions )23456789for iBeta = 1 : length ( betaList )betaDIR = [ DIR ’/ Delta = ’ num2str ( delta , ’ %.2 g ’) ...’/ Beta = ’ , num2str ( betaList ( iBeta ) , ’ %.2 g ’) ];if ~ isequal ( exist ( betaDIR , ’ dir ’) , 7) % 7 = directory .mkdir ( betaDIR ) ;endend1011121314numWorkers = feature ( ’ numCores ’) ;poolObj = parpool (2* numWorkers ) ;cleaner = onCleanup ( @ () delete ( poolObj ) ) ;parfor iBeta = 1 : length ( betaList )1516beta = betaList ( iBeta ) ;1718192021222324for iS = 1 : length ( sList )warning ( ’ off ’ , ’ MATLAB : odearguments : RelTolIncrease ’) ;s = sList ( iS ) ;betaDIR = [ DIR ’/ Delta = ’ num2str ( delta , ’ %.2 g ’) ’/ Beta = ’ ,num2str ( beta , ’ %.2 g ’) ];integrateTraj ( betaDIR , delta , beta , s , eqEpsVic , ODE , tEnd ,logOption , plotOptions ) ;warning ( ’ on ’, ’ MATLAB : odearguments : RelTolIncrease ’) ;end13225end2627figs2png ( delta , betaList , sList ) ;2829endЛистинг В.12:integrateTraj.m – функция, для интегрирования сепаратрисседловых состояний равновесия для заданных значений параметров , , .1function sol = integrateTraj ( DIR , delta , beta , s , eqEpsVic , ODE ,tEnd , logOption , plotOptions )23fileName = [ DIR , ’/s = ’ , num2str (s , ’ %.16 g ’) ];45warning ( ’ off ’ , ’ MATLAB : odearguments : RelTolIncrease ’) ;67S0 = [0 , 0 , 0]; S1 = [1 , 0 , 0]; % S2 = [ -1 , 0 , 0];8910% Values of parameters :alpha = delta * sqrt (1 - s ) ; lambda = s / sqrt (1 - s ) ;111213[V , D ] = eig ( J ( S0 , alpha , beta , lambda ) ) ;[ ~ , IX ] = sort ( real ( diag ( D ) ) , ’ descend ’) ;141516% % Eigenvectors of J in S0vU_S0 = V (: , IX (1) ) ’ / norm ( V (: , IX (1) ) ’) ;171819[ V1 , D1 ] = eig ( J ( S1 , alpha , beta , lambda ) ) ;[ ~ , IX1 ] = sort ( real ( diag ( D1 ) ) , ’ descend ’) ;202122if ~ isempty ( logOption )fprintf ( ’ %.12 g \ n ’ , s ) ;2324fprintf ( ’ parameters : alpha = %.16 g , beta = %.16 g , lambda = %.16g \ n ’ , alpha , beta , lambda ) ;2526272829fprintf ( ’ lambda_1 ( S0 ) = %s , lambda_2 ( S0 ) = %s , lambda_3 ( S0 ) = %s \ n ’ , D ( IX (1) , IX (1) ) , D ( IX (2) , IX (2) ) , D ( IX (3) , IX (3) ) ); %if real ( D( IX (1) , IX (1) ) ) > 0fprintf ( ’ Eq .
S0 is unstable \ n ’) ;else133fprintf ( ’ Eq . S0 is stable \ n ’) ;3031end3233343536373839fprintf ( ’ lambda_1 ( S1 ) = %s , lambda_2 ( S1 ) = %s , lambda_3 ( S1 ) = %s \ n ’ , D1 ( IX1 (1) , IX1 (1) ) , D1 ( IX1 (2) , IX1 (2) ) , D1 ( IX1 (3) , IX1 (3) ) ) ; %if real ( D1 ( IX1 (1) , IX1 (1) ) ) > 0fprintf ( ’ Eq . S1 and S2 are unstable \ n ’) ;elsefprintf ( ’ Eq . S1 and S2 are stable \n ’) ;endend404142434445464748495051initPointS0 = S0 + eqEpsVic .
S0 * sign ( vU_S0 (1) ) * vU_S0 ;sol . sepaS0 = feval ( ODE . solverName , @ (t , x ) lorenzLikeSyst (t , x ,alpha , beta , lambda ) , [0 , tEnd . trans ] , initPointS0 , ODE . options ) ;if ~ isempty ( plotOptions )plotLorenzLikeSepa ( fileName , sol . sepaS0 .x , sol . sepaS0 .y ’ ,plotOptions ) ;if ODE . limtLimSpan = linspace ( tEnd . trans - tEnd . lim , tEnd . trans ,100*tEnd . lim ) ;attrS0 = deval ( sol . sepaS0 , tLimSpan ) ;plotLimOptions = plotOptions ; plotLimOptions . color = false ;plotLorenzLikeSepa ( fileName , tLimSpan , attrS0 ’ ,plotLimOptions ) ;endend52535455if beta > s * ( delta + 2 / ( delta * (1 - s ) + s) )vU1_S1 = real ( V1 (: , IX1 (1) ) ) ’ / norm ( real ( V1 (: , IX1 (1) ) ) ) ;vU2_S1 = imag ( V1 (: , IX1 (1) ) ) ’ / norm ( imag ( V1 (: , IX1 (1) ) ) ) ;565758596061initPointS1 = S1 + eqEpsVic .