Диссертация (1145368), страница 25
Текст из файла (страница 25)
, 0; ...0, 0, 0, 1.; ...4.00 , 0 , -.40 , 0; ...0, -4. , 0, -.120];q = [0; 0; 12.6; 0.432];r = [ -6.0; -6.; -2.; -2.];% Parameter M in nonlinearity sat_low (z) = 1/2*( abs (z + M) - abs (z - M ))M = 5/57.3;% ODE solver parametersacc = 1e -13;RelTol = acc ; AbsTol = acc ; InitialStep = acc /10;options = odeset ( ' RelTol ', RelTol , ' AbsTol ', AbsTol , ' InitialStep ' , InitialStep ); % ,' NormControl ' ,'on '% Define harmonic linearization paramener k18326272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374I = eye ( size (A ));Wm = transpose (c )*( A -I*m*i )^( -1)* b;m_all = solve ( factor ( imag ( Wm )) , 'm ');% m1 = double (( m_all (3)))% k1 = - double (( real ( subs (Wm , 'm ', m1 )))^( -1))m2 = double (( m_all (5)));k2 = - double (( real ( subs (Wm , 'm ',m2 )))^( -1));% Define matrix P0P0_2 = P + k2 *q *r ';% P0_1 = P + k1 *q *r ';% Initial pointZ = [0.1 0.1 -0.2 -0.2];% Hidden oscillation localizationepsfig = figure ( ' Name ',' Hidden oscillation localization ');NStep = 10;for i = 1: NStepZsizes = size (Z ); lz = Zsizes (1);% Initial data for each new stepZ0 = Z(lz ,:);if ( abs ( Z0 )* abs ( Z0 ) ') > 1000returnend% Calculating eps valueeps = i/ NStep ;% Numerical integration of the system[T , Z] = ode45 ( @fSys_lin_in_eps ,[0 100] , Z0 , options , ...P0_2 , q , r , M , k2 , eps );end% Constructing plotsubplot (1 ,2 , mod (i +1 ,2)+1)plot (Z (: ,1) , Z (: ,2)); % ,'red ', ' blue ',' green 'grid on ; axis square ; hold on ;title ([ '\ epsilon = ' , num2str ( eps )]);xlabel ( ' x_1 '); ylabel ( ' x_2 '); zlabel ( ' x_3 ');if ( mod (i ,2)==0 )hold offdisp ( ' off ')h= figure ()endendÐåàëèçàöèÿ àëãîðèòìà âû÷èñëåíèÿ ëÿïóíîâñêîé ðàçìåðíîñòèàòòðàêòîðà Ýíîíà.Ðàññìîòðèì ìîäèôèêàöèþ àëãîðèòìà èç [188] äëÿ äèñêðåòíîé ñèñòåìûÝíîíà.Listing 3.4:1234567function out = henonMap ( x , a , b )out = zeros (6 ,1);out (1) = a + b*x (2) - x (1)^2;out (2) = x (1);out (3:6) = [ -2* x (1) , b; 1, 0];henonMap.m184İ = 0.2İ = 0.40.10.10.1−0.2−0.50−0.10x−0.2−0.50.510x020.1x20.220.2−0.1−0.10x−0.20.500.2x−0.2−0.50.41İ = 0.6İ = 0.7İ = 0.80.10.10.10.10x0.5−0.2−0.51−0.10x−0.2−0.50.5(c)0x0.5−0.2−0.51İ = 0.9ε = 0.7; 0.8İ = 1.00.20.10.120.2x0−0.10−0.10x0.5−0.2−0.510x0.51(e)0x1(d)ε = 0.5; 0.6−0.2−0.50−0.11x2−0.2−0.5−0.10xx−0.120.2x20.220.200.5ε = 0.3; 0.40.200x1(b)ε = 0.1; 0.2İ = 0.50−0.11(a)x2İ = 0.30.2xx2İ = 0.10.2ε = 0.9; 1.0Ðèñóíîê 3.25: Ïðîöåäóðà ëîêàëèçàöèè ñêðûòîãî êîëåáàíèÿ0.5185İ = 0.1İ = 1.00.10.1İȥ(ı)İȥ(ı)0.050.05ıı-0.05-0.05-0.1-0.1(a)(b)ε = 0.1ε = 1.0Ðèñóíîê 3.26: Ñåêòîð óñòîé÷èâîñòè íóëåâîãî ñîñòîÿíèÿ ðàâíîâåñèÿ èíåëèíåéíîñòü8endListing 3.5:12345678910111213141516171819202122232425262728293031323334henonAttrSE.mfunction henonAttrSEfunction out = henonMap ( x , a , b )out = zeros (2 ,1);endout (1) = a + b*x (2) - x (1)^2;out (2) = x (1);function out = J( x , a , b )out = [ -2* x (1) , b; 1, 0];enda = 1.4; b = 0.3;S1 = 1/2*(( b -1) + sqrt (( b -1)^2 + 4* a ));S2 = 1/2*(( b -1) - sqrt (( b -1)^2 + 4* a ));[V1 , D1 ] = eig (J ([ S1 , S1 ], a , b ));[~ , IX1 ] = sort ( diag ( D1 ), ' descend ');[V2 , D2 ] = eig (J ([ S2 , S2 ], a , b ));[~ , IX2 ] = sort ( diag ( D2 ), ' descend ');delta = 0.001;initPoint = [S2 , S2 ] + delta * V2 (: , IX2 (1)) ';nIter = 1000;sol = zeros ( nIter , 2);for i = 1 : nItersol (i , :) = henonMap ( initPoint , a , b ) ';initPoint = sol (i , :);endfigure (1)scatter ( sol (2: end , 1) , sol (2: end ,2) , ' filled ');18635363738394041endhold on ;plot (S2 , S2 , '.
', ' markersize ', 20 , ' Color ' , ' red ');grid on ;axis on ;xlabel ( 'x ' );ylabel ( 'y ' );Listing 3.6:123456789101112131415161718192021222324252627282930313233343536function [U , R , V] = productSVD ( initFactorization , nIterations )[~ , dimOde , nFactors ] = size ( initFactorization );A = zeros ( dimOde , dimOde , nFactors , nIterations );A (: , :, :, 1) = initFactorization ;Q = zeros ( dimOde , dimOde , nFactors +1);U = eye ( dimOde ); V = eye ( dimOde );R = zeros ( dimOde , dimOde , nFactors );for iIteration = 1 : nIterationsQ (: , :, nFactors + 1) = eye ( dimOde , dimOde );for jFactor = nFactors : -1 : 1C = A (: , :, jFactor , iIteration ) * Q (: , :, jFactor +1);[Q (: , : , jFactor ) , R (: , :, jFactor )] = qr (C );for kCoord = 1 : dimOdeif R( kCoord , kCoord , jFactor ) < 0R( kCoord , :, jFactor ) = -1 * R ( kCoord , : , jFactor );Q (: , kCoord , jFactor ) = -1 * Q (: , kCoord , jFactor );end ;end ;end ;if mod ( iIteration , 2) == 1U = U * Q (: , :, 1);elseV = V * Q (: , :, 1);endendendfor jFactor = 1 : nFactorsA (: , :, jFactor , iIteration + 1) = R (: , :, nFactors - jFactor +1) ';endListing 3.7:123456789101112131415161718productSVD.mcomputeLEsDiscret.mfunction LEs = computeLEsDiscret ( extMap , initPoint , ...dimOde = length ( initPoint );dimExtOde = dimOde * ( dimOde + 1);initCond = initPoint (:);X = zeros ( dimOde , dimOde , nFactors );for iFactor = 1 : nFactorsextMapSolution = extMap ( initCond );endX (: , :, iFactor ) = reshape (...extMapSolution (( dimOde + 1) : dimExtOde ), ...dimOde , dimOde );initCond = extMapSolution (1 : dimOde );18719202122232425262728[~ , R , ~] = productSVD (X , nSvdIterations );LEs = zeros (1 , dimOde );for jFactor = 1 : nFactorsLEs = LEs + log ( diag (R (: , :, jFactor )) ');end ;LEs = LEs / nFactors ;endListing 3.8:1234567891011121314151617181920lyapunovDim.mfunction LD = lyapunovDim ( LEs )LD = 0;nLEs = length ( LEs );sortedLEs = sort ( LEs , ' descend ');leSum = sortedLEs (1);if ( sortedLEs (1) > 0 )for i = 1 : nLEs -1if sortedLEs (i +1) ~= 0LD = i + leSum / abs ( sortedLEs (i +1) );leSum = leSum + sortedLEs (i +1);if leSum < 0break ;endendendendendListing 3.9:1234567891011121314151617main.mfunction maina = 1.4; b = 0.3;initPoint = [1;1];nFactors = 10000;nSvdIterations = 1;LEs = computeLEsDiscret (@( x) henonMap (x , a , b ), initPoint , nFactors , nSvdIterations );fprintf ( ' Lyapunov exponents : %6.4 f , %6.4 f , %6.4 f \n ', LEs );LD = lyapunovDim ( LEs );endfprintf ( ' Lyapunov dimension : %6.4 f\n ', LD );.