Диссертация (1144294), страница 16
Текст из файла (страница 16)
mat ’ , ’ traj ’) ;6768digits ( digitsOld ) ;6970end110Листинг Б.6:plotTraj2D.m – вспомогательная функция, рисующаятраекторию в двумерной плоскости.1function plotTraj2D ( traj , m1 , m2 ,b , coordInd2d )234syms t positivesyms X_1 X_2 X_3 X_4 real56[ numArcs , ~] = size ( traj ) ;78figure (1) ; hold on ;910color = { ’ red ’, ’ blue ’} ;111213for iTr = 1 : numArcssymSolOde = str2func ([ ’ symSolOde ’ , int2str ( traj { iTr , 1 }) ]) ;14currSol = symSolOde (m1 , m2 , b) ;currIC = traj { iTr , 2 } ;sol_x1 = subs ( currSol ( coordInd2d { 1 } ) , { X_1 , X_2 , X_3 , X_4 } , ...{ currIC (1) , currIC (2) , currIC (3) , currIC (4) } ) ;sol_x2 = subs ( currSol ( coordInd2d { 2 } ) , { X_1 , X_2 , X_3 , X_4 } , ...{ currIC (1) , currIC (2) , currIC (3) , currIC (4) } ) ;1516171819202122232425262728plot3 ( currIC ( coordInd2d {1 } ) , currIC ( coordInd2d { 2 } ) , ’.
’ , ’markersize ’ , 15 , ’ Color ’ , ’ black ’) ;fplot ( sol_x1 , sol_x2 , [0 , double ( traj { iTr , 3 } ) ] , ’ Color ’ , color {traj { iTr , 1} } ) ;endgrid on ; axis on ;xlabel ([ ’x ’ int2str ( coordInd2d { 1 } )] , ’ FontName ’ , ’ Times New Roman ’ ,’ FontAngle ’ , ’ italic ’, ’ FontSize ’ , 21) ;ylabel ([ ’x ’ int2str ( coordInd2d { 2 } )] , ’ FontName ’ , ’ Times New Roman ’ ,’ FontAngle ’ , ’ italic ’, ’ FontSize ’ , 21) ;hold off ;2930endЛистинг Б.7:plotTraj3D.m – вспомогательная функция, рисующаятраекторию в трехмерной гиперплоскости.1111function plotTraj3D ( traj , m1 , m2 , b , coordInd3d )234syms t positivesyms X_1 X_2 X_3 X_4 real56[ numArcs , ~] = size ( traj ) ;78figure (1) ; hold on ;910trajColors = { ’ red ’ , ’ blue ’} ;1112131415161718192021for iTr = 1 : numArcssymSolOde = str2func ([ ’ symSolOde ’ , int2str ( traj { iTr , 1 }) ]) ;currSol = symSolOde (m1 , m2 , b) ;currIC = traj { iTr , 2 } ;sol_x1 = subs ( currSol ( coordInd3d { 1 } ) , { X_1 , X_2 , X_3 , X_4 } , ...{ currIC (1) , currIC (2) , currIC (3) , currIC (4) } ) ;sol_x2 = subs ( currSol ( coordInd3d { 2 } ) , { X_1 , X_2 , X_3 , X_4 } , ...{ currIC (1) , currIC (2) , currIC (3) , currIC (4) } ) ;sol_x3 = subs ( currSol ( coordInd3d { 3 } ) , { X_1 , X_2 , X_3 , X_4 } , ...{ currIC (1) , currIC (2) , currIC (3) , currIC (4) } ) ;222324% Plot switching pointplot3 ( currIC ( coordInd3d { 1 } ) , currIC ( coordInd3d { 2 }) , currIC (coordInd3d { 3 } ) , ’.
’ , ’ markersize ’ , 15 , ’ Color ’, ’ black ’) ;25262728% Plot the arc of trajectory of the regimefplot3 ( sol_x1 , sol_x2 , sol_x3 , [0 , double ( traj { iTr , 3} ) ] , ’Color ’ , trajColors { traj { iTr , 1 } } );end29303132333435grid on ; axis on ;xlabel ([ ’x ’ int2str ( coordInd3d { 1 } )] , ’ FontName ’ , ’ Times New Roman ’ ,’ FontAngle ’ , ’ italic ’, ’ FontSize ’ , 21) ;ylabel ([ ’x ’ int2str ( coordInd3d { 2 } )] , ’ FontName ’ , ’ Times New Roman ’ ,’ FontAngle ’ , ’ italic ’, ’ FontSize ’ , 21) ;zlabel ([ ’x ’ int2str ( coordInd3d { 3 } )] , ’ FontName ’ , ’ Times New Roman ’ ,’ FontAngle ’ , ’ italic ’, ’ FontSize ’ , 21) ;hold off ;11236endЛистинг Б.8:findPeriodicExact.m – функция, для аналитико-численнойлокализации периодического решения методом точечных отображений.1function [ solX1 , solX2 , solX4 , solT1 , solT2 ] = findPeriodicExact ( m1 , m2 ,beta , INIT_GUESS )2345syms t t1 t2 positivesyms X_1 X_2 X_3 realsyms X_4 positive678solOde1 = symSolOde1 ( m1 , m2 , beta ) ;solOde2 = symSolOde2 ( m1 , m2 , beta ) ;910111213x1_Ode1x2_Ode1x3_Ode1x4_Ode1====subs ( solOde1 (1) ,subs ( solOde1 (2) ,subs ( solOde1 (3) ,subs ( solOde1 (4) ,{t ,{t ,{t ,{t ,X_3 } ,X_3 } ,X_3 } ,X_3 } ,{ -t1 ,{ -t1 ,{ -t1 ,{ -t1 ,x1_Ode2x2_Ode2x3_Ode2x4_Ode2====subs ( solOde2 (1) ,subs ( solOde2 (2) ,subs ( solOde2 (3) ,subs ( solOde2 (4) ,{t ,{t ,{t ,{t ,X_3 } ,X_3 } ,X_3 } ,X_3 } ,{ t2 ,{ t2 ,{ t2 ,{ t2 ,0});0});0});0});14151617180});0});0});0});192021222324assume ( X_4 > 0) ; assume ( t1 > 0) ; assume ( t2 > 0) ;[ solX1 , solX2 , solX4 , solT1 , solT2 ] = vpasolve ( ...[ x1_Ode2 - x1_Ode1 == 0 , x2_Ode2 - x2_Ode1 == 0 , ...x3_Ode2 - x3_Ode1 == 0 , x3_Ode2 == 0 , x4_Ode2 - x4_Ode1 == 0] , [X_1 , X_2 , X_4 , t1 , t2 ] , INIT_GUESS ) ;endЛистинг Б.9:fittsMain.m – функция, для запуска аналитико-численногоалгоритма поиска периодической траектории системы (1.21) параметрами1 = 0.9, 2 = 1.1, = 0.03.1function fittsMain23vpaPresision = 32;45m1 = vpa ( ’ 0.9 ’ , vpaPresision );11367m2 = vpa ( ’ 1.1 ’ , vpaPresision );beta = vpa ( ’ 0.03 ’ , vpaPresision );8910x0 = [ vpa ( ’ 10 ’ , vpaPresision ); vpa ( ’ 10 ’ , vpaPresision ) ;vpa ( ’ 10 ’ , vpaPresision ) ; vpa ( ’10 ’ , vpaPresision ) ];1112tEnd = vpa ( ’ 500 ’ , vpaPresision ) ;1314traj = integrateTrajectory (m1 , m2 , beta ,x0 , tEnd , vpaPresision );151617plotTraj3D ( traj , m1 , m2 , b , { 1 ,2 ,3 }) ;plotProj ( traj , m1 , m2 , b , { 1 ,2 } ) ;181920INIT_GUESS = [ traj { end -1 ,2 } (1) , traj { end -1 ,2 } (2) , ...traj { end -1 ,2 } (4) , traj { end -2 ,3 } (1) , traj { end -1 ,3 } (1) ];2122[ solX1 , solX2 , solX4 , solT1 , solT2 ] = findPeriodicExact ( m1 , m2 , beta ,INIT_GUESS ) ;232425disp ( solX1 ) ; disp ( solX2 ) ; disp ( solX4 ) ;disp ( solT1 ) ; disp ( solT2 ) ;2627endБ.2Локализация хаотического решения в системе Фиттса с разрывной правой частью методом продолженияпо параметруВоспользуемся для интегрирования решений системы (1.21) с разрывнойправой частью специальным вычислительной процедурой filippov() для численного моделирования решений по Филиппову [93].Листинг Б.10:fittsVectorFields.m – функция, задающая правые частисистемы (1.21) в областях Σ± .1function [F1 , F2 ,H , dH ,h , dir ] = fittsVectorFields (t ,y , params , str )23% Parameters :114456m1 = params (1) ;m2 = params (2) ;beta = params (3) ;7891011a3a2a1a0====4* beta ;m1 ^2+ m2 ^2+6* beta ^2;(2*( m1 ^2+ beta ^2) ) * beta +2* beta *( m2 ^2+ beta ^2) ;( m1 ^2+ beta ^2) *( m2 ^2+ beta ^2) ;12131415% Vector field in region 1 - H ( x ) > 0F1 = [ y (2) ; y (3) ; y (4) ;- a3 * y (4) - a2 * y (3) - a1 * y (2) - a0 * y (1) + 1];16171819% Vector field in region 2 - H ( x ) < 0F2 = [ y (2) ; y (3) ; y (4) ;- a3 * y (4) - a2 * y (3) - a1 * y (2) - a0 * y (1) - 1];202122% Switching ManifoldH = -y (3) ;232425% A vector normal to the switching manifolddH = -[0 ,0 ,1 ,0];262728% Poincare sectionh = 1;29303132% Location directriondir = 1;endЛистинг Б.11:fittsJacobians.m – функция, задающая линейные частисистемы (1.21) в областях Σ± .1function [J1 , J2 , d2H ] = fittsJacobians (t ,y , params , str )23456% Parametersm1 = params (1) ;m2 = params (2) ;beta = params (3) ;78a3 = 4* beta ;11591011a2 = m1 ^2+ m2 ^2+6* beta ^2;a1 = (2*( m1 ^2+ beta ^2) ) * beta +2* beta *( m2 ^2+ beta ^2) ;a0 = ( m1 ^2+ beta ^2) *( m2 ^2+ beta ^2) ;121314151617% Jacobian in region S1 ( H( x ) > 0)J1 = [0 , 1 , 0, 0;0 , 0 , 1 , 0;0 , 0 , 0 , 1;-a0 , -a1 , -a2 , - a3 ];181920% Jacobian in region S2 ( H( x ) < 0)J2 = J1 ;21222324% grad ( grad ( H ) ) A vector normal to the discontinuity surfaced2H = zeros ( size ( J1 ) ) ;endЛистинг Б.12:fittsContinuation.m – функция, для локализацияхаотического решения в системе Фиттса с разрывной правой частью методомпродолжения по параметру.1function fittsContinuation234% ODE solversolver = ’ ode45 ’;567% ODE - solver optionsopts = odeset ( ’ RelTol ’ ,1e -8 , ’ AbsTol ’ ,1e -8 , ’ MaxStep ’ ,0.1) ;8910% Name of the file with the two vectorfieldsvfields = ’ fittsVectorFields ’;111213% Name of the file with the two Jacobianssjacobians = ’ fittsJacobians ’;141516% Name of the file with the two Jacobiansspfunction = ’ ’;171819% Filippov parameterC = 1;11620212223% Parametersm1 = 0.9; m2 = 1.1;beta0 = 0.03; betaStep = 0.0175; nSteps = 4;242526% Actual integration timeT = 2000; tSpan = [0 , T ];2728293031% Initial conditions of the states obatined analyticallyx0 = [ -0.62520516260693109534342362490723 , ...-3.7324097072650610465825278562594 , 0 , ...3.4754169728697120793989274111636];323334for iStep = 1 : nStepsbeta = beta0 + ( iStep -1) * betaStep ;3536params = [ m1 , m2 , beta ];37383940% Output is the time , states and events as in Matlab ’s standardoutput[ ~ ,x , ~ ,~ ,~ ,~ ] = filippov ( vfields , jacobians , pfunction , solver ,tSpan , x0 , params ,C , opts ) ;x0 = x ( end , :) ;41424344454647hFig1 = figure ; hold on ;plot3 ( x (: ,1) , x (: ,2) , x (: ,3) , ’ Color ’ , [0 , 0 , 1]) ;axis on ; grid on ;xlabel ( ’ x_1 ’) , ylabel ( ’ x_2 ’) , zlabel ( ’ x_3 ’)view ([42 ,32]) ; hold off ; box off ;savefig ( hFig1 ,[ ’ fitts_cont_ ’ int2str ( iStep -1) ’ _x1x2x3 .