programming to controller design for discrete input hybrid nonlinear
A.1 The Mathematica package
(∗ : : Package : : ∗) (∗ : : T i t l e : : ∗) (∗Maximal Lyapunov∗)
(∗ : : Section : : ∗)
(∗BeginPackage [" MaximalLyapunov ` " ]∗)
BeginPackage [ "MaximalLyapunov ` " ]
Grad : : usage = " Gives the g r a d i e n t vector o f i t s f i r s t argument , r e s p e c t to the second one . "
NormalizeMatrix : : usage = " Normalizes the rows o f the matrix given as parameter . "
SquareMatrixQ : : usage = "Checks i f a matrix i s square matrix . "
IsLocalAsymptoticalStable : : usage = "Checks i f the system d e s c r i b e d by the matrix given as parameter i s l o c a l l y asymptotical s t a b l e . "
SymmetricalMatrix : : usage = " Generates a symmetrical matrix . "
I f [ $VersionNumber<=7,
LyapunovSolve : : usage = " Solves the Lyapunov−equation . "
MultipleLimit : : usage = " C a l c u l a t e s the limes along m u l t i p l e v a r i a b l e s . "] ;
PiecewiseExtract : : usage = " Extracts the function−parts from e x p r e s s i o n s o f the form o f Piecewise [ \ [ E l l i p s i s ] ] "
PiecewiseConditionExtract : : usage = " Extracts the condition−parts from e x p r e s s i o n s o f the form o f Piecewise [ \ [ E l l i p s i s ] ] "
PiecewiseCombine : : usage = "Combines p i e c e w i s e e x p r e s s i o n s i n t o one . "
InitMaximalLyapunov : : usage = " I n i t i a l i z e s the Vanelli−Vidyasagar algorithm . "
IterateMaximalLyapunov : : usage = " Executes one i t e r a t i o n o f the Vanelli− Vidyasagar algorithm . Parameters can be any o f what can be used in
f u n c t i o n NMinimize furthermore Chop can be e i t h e r True or False . "
MaximalLyapunovFunction : : usage = " Represents the maximal Lyapunov f u n c t i o n . "
GradientOfMaximalLyapunovFunction : : usage = " Representes the g r a d i e n t o f the maximal Lyapunov f u n c t i o n . "
LyapunovQs : : usage = " Gives the sum o f Q' s in the estimated maximal Lyapunov f u n c t i o n . "
L i s t O f V a r i a b l e s : : usage = " Contains the s t a t e v a r i a b l e s f o r the i n i t i a l i z e d system . "
PlotDOA : : usage = " Plots the DOA o f the system . "
EstimateIntervalOfCStar : : usage = " This f u n c t i o n e s t i m a t e s an i n t e r v a l which \ ! \ ( \∗SuperscriptBox [ \ (C\) , \(∗\) ] \ ) l i e s in . "
ConvergesToTheOriginQ : : usage = "Checks i f the provided s o l u t i o n converges to the o r i g i n . "
ScanDomainOfStability : : usage = " Scans the DOA o f the by brute f o r c e . "
(∗ : : Section : : ∗) (∗Begin of ` Private `∗)
Begin [ " ` Private ` " ]
(∗ Mathematica Raw Program ∗)
(∗ : : Subsection : : Closed : : ∗) (∗Vector and Matrix o p e r a t i o n s∗)
SyntaxInformation [ Grad ] = {" ArgumentsPattern "−>{_, _} } ; Grad [ f_ , vars_List ] :=
D[ f , { vars } ]
SyntaxInformation [ NormalizeMatrix ] = {" ArgumentsPattern "−>{_} } ; NormalizeMatrix [x_?MatrixQ ] :=
Map[#/Norm[#]& , x ]
SyntaxInformation [ IsLocalAsymptoticalStable ] = {" ArgumentsPattern "−>{_
} } ;
IsLocalAsymptoticalStable [x_?MatrixQ ] :=
And@@((#<0&)/@Re[ Eigenvalues [ x ] ] )
SyntaxInformation [ SymmetricalMatrix ] = {" ArgumentsPattern "−>{_, _} } ; SymmetricalMatrix [ n_Integer , h_Symbol ] :=
Array [ h , {n , n } ] / . m_[ r_ , c_ ] :> m[ c , r ] / ; r > c
SyntaxInformation [ SquareMatrixQ ] = {" ArgumentsPattern "−>{_} } ; SquareMatrixQ [ mat_ ] :=
MatrixQ [ mat ] && Equal @@ Dimensions [ mat ]
I f [ $VersionNumber<=7,
(∗ " CompatibilityMode " can be e i t h e r " Matlab " or "Mathematica" ∗) Options [ LyapunovSolve ] = {" SymmetricForce " −> False , "
CompatibilityMode " −> "Mathematica" } ;
SyntaxInformation [ LyapunovSolve ] = {" ArgumentsPattern "−>{_, _, OptionsPattern [ ] } } ;
LyapunovSolve [A_? SquareMatrixQ , Q_? SquareMatrixQ , OptionsPattern [ ] ] Module [ {X, m} ,:=
m = I f [ OptionValue [ " SymmetricForce " ] ,
SymmetricalMatrix [ Dimensions [A ] [ [ 1 ] ] , x ] , Array [ x,{#,#}&@Dimensions [A ] [ [ 1 ] ] ]
X = (#/. Solve [A.#+#.Transpose [A]==] ; −Q, Union [ Flatten [ # ] ] ] [ [ 1 ] ] ) &[m] ;
I f [ OptionValue [ " CompatibilityMode " ] == "Matlab" , Return [X] ,
Return[−X]
] ] ]
SyntaxInformation [ MultipleLimit ] = {" ArgumentsPattern "−>{_, _} } ; MultipleLimit [ f_ , r u l e s _ L i s t ] :=
Module [ {numer , tog } , tog = Together [ f ] ;
numer = Numerator [ tog ] / . r u l e s ; tog = numer/Denominator [ tog ] ; Return [ tog / . r u l e s ] ;
]
(∗ : : Subsection : : Closed : : ∗)
(∗Operations on Piecewise f u n c t i o n s∗)
PiecewiseExtract [ HoldPattern [ Piecewise [ l :{ {_, _} . . } , last___ ] ] ] :=
Flatten [ { l [ [ All , 1 ] ] , l a s t } ] PiecewiseExtract [x_] = {x}
PiecewiseConditionExtract [ HoldPattern [ Piecewise [ l :{ {_,_} . . } , last___ ] ] ] Flatten [ { l [ [ All , 2 ] ] , True } ]:=
PiecewiseConditionExtract [_] :=
{True}
SyntaxInformation [ PiecewiseCombine ] = {" ArgumentsPattern "−>{_, _} } ; PiecewiseCombine [ body_List , condition_List ] :=
Piecewise [ Inner [ List , body , condition , List ] ]
(∗ : : Subsection : : Closed : : ∗)
(∗S t r u c t u r a l o p e r a t i o n s and f u n c t i o n s∗)
GetTag : : usage = "Get the name o f the v a r i a b l e s in a polyonimial e x p r e s s i o n . "
SyntaxInformation [ GetTag ] = {" ArgumentsPattern "−>{_, _} } ; GetTag [ expr_ , coeff_Symbol ] :=
Cases [ Variables [ PiecewiseExtract @ PiecewiseExpand [ expr ] ] , _coeff | c o e f f ]
SyntaxInformation [ M u l t i p l i e r L i s t ] = {" ArgumentsPattern "−>{_, _, _} } ; M u l t i p l i e r L i s t [ expr_ , t_List , v_List ] :=
Coefficient [ expr , t ] / .Map[ Rule[# , 0]& , v ]
SyntaxInformation [ GetTermsGreaterThanOrEqual ] = {" ArgumentsPattern "−>{_
, _, _} } ;
GetTermsGreaterThanOrEqual [ expr_ , n_Integer , x_Symbol ] :=
Pick [ expr ,
Map[#>=n&, Plus@@(Exponent [ List@@expr , #]&/@GetTag [ expr , x ] ) ] ]
SyntaxInformation [ VariableOrderList ] = {" ArgumentsPattern "−>{_, _} } ; VariableOrderList [ i_Integer , n_Integer ] :=
Reverse [ Select [ Tuples [ Range[ 0 , i ] , n ] , Plus@@#===i &]]
GenerateMonomialList : : usage = " Creates the l i s t o f m u l t i p l i c a n d s in a monomial o f order i and with v a r i a b l e s o f number n . "
SyntaxInformation [ GenerateMonomialList ] = {" ArgumentsPattern " −> {_, _ GenerateMonomialList [ i_Integer , n_Integer ] [ x_Symbol ] :=}}
Map[ Flatten [MapThread[ Table[#1 , {#2}]&, {Array [ x , n ] , #}]]& , VariableOrderList [ i , n ] ]
SyntaxInformation [ Monomials ] = {" ArgumentsPattern " −> {_, _}}
Monomials [ i_Integer , n_Integer ] [ x_Symbol ] :=
Times@@@GenerateMonomialList [ i , n ] [ x ]
(∗ : : Subsection : : ∗)
(∗Approximation of maximal Lyapunov f u n c t i o n∗)
SyntaxInformation [ Fi ] = {" ArgumentsPattern " −> {_, _}}
Fi [ i_Integer ? Positive , n_Integer ? Positive ] [ x_Symbol , a_Symbol ] :=
Plus@@MapThread [ Times , {#, Table [ a [ k , i ] , {k , Length[#]}]}&
@Monomials [ i , n ] [ x ] ] Fi [ _Integer , _Integer ] [__] = 0 ;
V : : usage = " Represents the maximal Lyapunov f u n c t i o n . "
SyntaxInformation [V] = {" ArgumentsPattern " −> {_, _}}
V[ i_Integer , varnum_Integer ] [ x_Symbol ] :=
Sum[R[ j , varnum ] [ x ] , { j , 2 , i }]/(1+Sum[Q[ j , varnum ] [ x ] , { j , 1 , −2+i } ] )
SyntaxInformation [R] = {" ArgumentsPattern " −> {_, _}}
R[ n_Integer ? Positive , varnum_Integer ? Positive ] [ x_Symbol ] :=
Fi [ n , varnum ] [ x , Global ` a ] R[ _Integer , _Integer ] [ _Symbol ] = 0 ;
SyntaxInformation [Q] = {" ArgumentsPattern " −> {_, _}}
Q[ n_Integer ? Positive , varnum_Integer ? Positive ] [ x_Symbol ] :=
Fi [ n , varnum ] [ x , Global ` b ] Q[ _Integer , _Integer ] [ _Symbol ] = 0 ;
Phi : : usage = "Option ReplaceFunction can be ReplaceAll or MultipleLimit
"
Options [ Phi ] = { ReplaceFunction−>ReplaceAll}
SyntaxInformation [ Phi ] = {" ArgumentsPattern "−>{_, _, OptionsPattern [ ] } } ;
Phi [ k_Integer , centralODE_ , OptionsPattern [ ] ] :=
Module [ {nomin , denom , arg1 , arg2 } ,
nomin = Map[D[ centralODE , Sequence@@#1]&, GenerateMonomialList [ k , nov ] [ x ] ] ;
(∗ nomin = nomin /. Indeterminate−>OptionValue ["
ReplaceIndeterminate " ] ;∗)
denom = Times@@@Factorial [ VariableOrderList [ k , nov ] ] ; arg1 = Transpose [ nomin/denom ] ;
arg2 = Map[ Rule[# , 0]& , xvar [ ] ] ;
OptionValue [ ReplaceFunction ] [ arg1 , arg2 ] ]
F [ k_Integer?(#<=0&) , centralODE_ ] [ x_Symbol ] :=
Table [ 0 , {nov } ]
F [ k_Integer , centralODE_ ] [ x_Symbol ] :=
Phi [ k , centralODE ] . Monomials [ k , nov ] [ x ]
VnDot [ n_Integer , centralODE_ ] [ x_Symbol ] :=
Grad [V[ n , nov ] [ x ] , xvar [ ] ] . centralODE
U[ n_Integer ] [ x_Symbol ] :=
(1+Sum[Q[ i , nov ] [ x ] , { i , 1 , −2+n } ] ) ^2
GreatTerms [ n_Integer , centralODE_ ] [ x_Symbol ] :=
GetTermsGreaterThanOrEqual [Expand[ Cancel [ VnDot [ n , centralODE ] [ x ]U[ n ] [ x ] ] ] , n+1, x ]
GreatTermsHybrid [ n_Integer , centralODE_ ] [ x_Symbol ] :=
Expand[
PiecewiseExpand [
VnDot [ n , centralODE ] [ x ]∗U[ n ] [ x]−Grad [R[ 2 , nov ] [ x ] , xvar [ ] ] . F[ 1 , centralODE ] [ x]−
Sum[Grad [R[ 2 , nov ] [ x ] , xvar [ ] ] . F [ k−1, centralODE ] [ x]+
Sum[( Grad [R[ j , nov ] [ x ] , xvar [ ] ] +
Sum[Q[ i , nov ] [ x ]∗Grad [R[ j−i , nov ] [ x ] , xvar [ ] ]−Grad [Q [ i , nov ] [ x ] , xvar [ ] ]∗R[ j−i , nov ] [ x ] ,
{ i , 1 , j−2}]) . F [ k−j +1, centralODE ] [ x ] , { j , 3 , k } ] ,
{k , 3 , n } ] ] ]
SyntaxInformation [ e0 ] = {" ArgumentsPattern "−>{_} } ; e0 [ expr_ ] :=
Plus@@(#^2&/@List@@expr )
eNonHybrid [ n_Integer , centralODE_ ] [ x_Symbol ] :=
e0 [ GreatTerms [ n , centralODE ] [ x ] ]
eHybrid [ n_Integer , centralODE_ ] [ x_Symbol ] :=
Module [ { funpart , condpart , expr } ,
expr = GreatTermsHybrid [ n , centralODE ] [ x ] ; funpart = e0 /@Expand [ PiecewiseExtract [ expr ] ] ; condpart = PiecewiseConditionExtract [ expr ] ; PiecewiseCombine [ funpart , condpart ]
]
e q 4 4 l e f t [ 2 , centralODE_ ] [ x_Symbol ] :=
Grad [R[ 2 , nov ] [ x ] , xvar [ ] ] . F[ 1 , centralODE ] [ x ] e q 4 4 l e f t [ k_Integer?(#>=3&) , centralODE_ ] [ x_Symbol ] :=
Expand[Sum[ Grad [R[ i , nov ] [ x ] , xvar [ ] ] . F [ k+1−i , centralODE ] [ x ] , { i , 2 , k}]+
Sum[ (Q[ i , nov ] [ x ]∗Grad [R[ j , nov ] [ x ] , xvar [ ] ]−
Grad [Q[ i , nov ] [ x ] , xvar [ ] ]∗R[ j , nov ] [ x ] ) . F [ k+1−i−j , centralODE ] [ x ] ,
{ i , 1 , k−2}, { j , 2 , k−1}]]
e q 4 4 r i g h t [ 2 ] [ x_Symbol ] :=
−xvar [ ] . \ [ GothicCapitalQ ] . xvar [ ] e q 4 4 r i g h t [ k_Integer?(#>=3&) ] [ x_Symbol ] :=
Expand[−xvar [ ] . \ [ GothicCapitalQ ] . xvar [ ] ( 2Q[ k−2, nov ] [ x]+
Sum[Q[ i , nov ] [ x ]∗Q[ k−2−i , nov ] [ x ] , { i , 1 , k−3}]) ]
SyntaxInformation [ myColorPrint ] = {" ArgumentsPattern "−>{_, _} } ; myColorPrint [ string_String , value_Symbol ] :=
CellPrint [
TextCell [Row[ { s t r i n g , " \ [ I n v i s i b l e S p a c e ] " ,
S t y l e [ ToString@value , FontColor −> I f [ value , Blue , ] ] } ] , " Print " ,Red Editable −> False ] ]
Options [ InitMaximalLyapunov ] = {Verbose −> 1}
SyntaxInformation [ InitMaximalLyapunov ] = {" ArgumentsPattern " −> {_, _, _, _, OptionsPattern [ ] } } ;
InitMaximalLyapunov [ odeqr_List , equp_List , zz_Symbol , xx_Symbol , OptionsPattern [ ] ] :=
Module [ {centralODE , isHybrid , localAsStable , f i r s t P h i , P} , Clear [ Global ` a ] ;
Clear [ Global ` b ] ; Clear [ c o e f f s ] ; z = zz ;
x = xx ;
(∗ Check i f equp r e a l l y i s an e q u i l i b r i u m poi nt ∗) I f [ OptionValue [ Verbose]>=1,
Print [ "Check i f the given e q u i l i b r i u m i s a r e a l one : " , odeqr / . z [ i_ ]: > equp [ [ i ] ] ]
] ;
(∗ After the transformationx_i=z_i−z_i^∗we g e t : ∗)
centralODE = Together [Expand[ odeqr / . z [ i_Integer ]: > x [ i ]+equp [ [ i ] ] ] ] ;
I f [ OptionValue [ Verbose]>=1,
Print [ "The system a f t e r r e c e n t r a l i z a t i o n i s : \ n" , TableForm [ centralODE ] ]
] ;
(∗ Check i f the transformed system has hybrid behaviour ∗) isHybrid = ! FreeQ [ centralODE , Piecewise ] ;
I f [ OptionValue [ Verbose]>=1,
Print [ "The given system has hybrid behaviour : " , isHybrid ] eTag = I f [ isHybrid ,] ;
eHybrid , eNonHybrid
] ;
(∗ Get the number of v a r i a b l e s in the system being considered
∗)
(∗ nov = Length [ Union [ F l a t t e n [ Table [ GetTag [ # [ [ i ] ] , x ] , { i , Length [#]}]&[ centralODE ] ] ] ] ; ∗)
nov = Length [ equp ] ;
I f [ OptionValue [ Verbose]>=1,
Print [ "The number o f v a r i a b l e s in the system : " , nov ] ] ;
(∗ I n i t i a l i z e s v a r i a b l e s seen o u t s i d e of InitMaximalLyapunov ∗)
\ [ GothicCapitalQ ] = IdentityMatrix [ nov ] ; xvar [ ] = Table [ x [ i ] , { i , nov } ] ;
xvarunit = (#−>1)&/@xvar [ ] ; itnum = 2 ;
(∗ Check i f the system i s l o c a l l y a s y m p t o t i c a l l y s t a b l e ∗) f i r s t P h i = Phi [ 1 , centralODE ] ;
l o c a l A s S t a b l e = IsLocalAsymptoticalStable [ f i r s t P h i ] ; I f [ OptionValue [ Verbose]>=1,
myColorPrint [ "The given system i s l o c a l l y a s y m p t o t i c a l l y s t a b l e : " , l o c a l A s S t a b l e ]
I f [ ! localAsStable ,] ; Return [ { centralODE } ] (∗ ] ;
I f [ nov==2,
Module [ {x , y , eigvec , normEigVec } , e i g v e c = Eigenvectors [ f i r s t P h i ] ; normEigVec = Norm /@ e i g v e c ; Print [
Show [{
StreamPlot [ Evaluate [ f i r s t P h i .{ x , y } ] , { x , −1, 1} , {y , −1, 1}] ,
Graphics [{ Red , Arrow [1/2 Normalize /@
{{0 , 0} , #}]}& /@ e i g v e c ] ,
Graphics [{ Green , Arrow [1/2 {{0 , 0} , #}
/ Max[ normEigVec ]]}& /@ e i g v e c ] ] }]
] ; ]
∗)
P = −LyapunovSolve [ Transpose [ f i r s t P h i ] , \ [ GothicCapitalQ ] ] ; I f [ OptionValue [ Verbose]>=1,
Print [ "The eigensystem i s : " , Eigensystem [ f i r s t P h i ] / / MatrixForm ] ;
Print [ " Matrix P i s : " , P//MatrixForm ] ; c o e f f s =] ;
MapThread[
Rule , { M u l t i p l i e r L i s t [V[ itnum , nov ] [ x ] ,# , xvar [ ] ] ,
M u l t i p l i e r L i s t [Expand[ xvar [ ] . P. xvar [ ] ] , # , xvar [ ] ] } &@
Monomials [ itnum , nov ] [ x ] ] ; (∗ Quick check ∗)
I f [ OptionValue [ Verbose]>=1,
Print [ "Quick check 1 ; should be zero : " ,
Transpose [ f i r s t P h i ] . P+P. f i r s t P h i +\[ GothicCapitalQ ] / / MatrixForm ] ;
Print [ "Quick check 2 ; should be zero : " ,
Expand[ e q 4 4 l e f t [ 2 , centralODE ] [ x]−e q 4 4 r i g h t [ 2 ] [ x ] / . c o e f f s ] ] ;
] ;
{centralODE , itnum , c o e f f s } ]
(∗ As option "Chop−>True | False " and any option accepted by NMinimize can be given ∗)
Options [ IterateMaximalLyapunov ] = Union [ Options [ NMinimize ] , Options [ Minimize ] , {Chop −> False , Verbose −> 0 , MinimizerApproach −> "
Numeric" } ]
SyntaxInformation [ IterateMaximalLyapunov ] = {" ArgumentsPattern " −> {_, OptionsPattern [ ] } } ;
IterateMaximalLyapunov [ centralODE_ , opts : OptionsPattern [ ] ] :=
Module [ {y , l s i d e , bn , AnEb, mincrit , minsol } , (∗ i n c r e a s e the i t e r a t i o n counter ∗)
itnum++;
(∗ v a r i a b l e s ∗)
y = Union [ GetTag [Q[ itnum−2, nov ] [ x ] , Global ` b ] , GetTag [R[ itnum , nov ] [ x ] , Global ` a ] ] ;
l s i d e = M u l t i p l i e r L i s t [ e q 4 4 l e f t [ itnum , centralODE ] [ x ] / . c o e f f s , Monomials [ itnum , nov ] [ x ] ,
xvar [ ] ] ;
bn = M u l t i p l i e r L i s t [ e q 4 4 r i g h t [ itnum ] [ x ] / . c o e f f s , Monomials [ itnum , nov ] [ x ] ,
xvar [ ] ] ;
(∗ minimization c o n s t r a i n t s ∗)
AnEb = MapThread[ Equal , { l s i d e , bn } ] ; I f [ OptionValue [ Verbose]>=1,
Print [ "Number o f i t e r a t i o n s : " , itnum ] ; Print [ " Equations : " , AnEb ] ;
] ;
(∗ f u n c t i o n to be minimized ∗)
mincrit = eTag [ itnum , centralODE ] [ x ] / . c o e f f s / . xvarunit ; I f [ OptionValue [ Verbose]>=1,
Print [ " C r i t e r i a : " , mincrit ] ; Print [ " V a r i a b le s : " , y ] ; ] ;
(∗ do the minimization ∗)
I f [ OptionValue [ MinimizerApproach ] == "Numeric" ,
minsol = NMinimize[ { mincrit , AnEb} , y ,
Sequence@@FilterRules [ { opts } , Options@NMinimize ] ] , minsol = Minimize [ { mincrit , AnEb} , y ,
Sequence@@FilterRules [ { opts } , Options@Minimize ] ] I f [ OptionValue [ Verbose]>=1,] ;
Print [ " S o l u t i o n : " , minsol ] ; c o e f f s = Union [ c o e f f s , minsol [ [ 2 ] ] ] ;] ; I f [ OptionValue [Chop] ,
c o e f f s = Chop@coeffs ] ;
MaximalLyapunovFunction = V[ itnum , nov ] [ x ] / . c o e f f s ;
GradientOfMaximalLyapunovFunction = VnDot [ itnum , centralODE ] [ x ] / . c o e f f s ;
Return [ { itnum , minsol [ [ 1 ] ] } ] ; ]
(∗ : : Subsection : : Closed : : ∗) (∗R e s u l t s of the algorithm∗) LyapunovQs :=
Sum[Q[ n , nov ] [ x ] , {n , 1 , itnum−2}]/. c o e f f s
L i s t O f V a r i a b l e s :=
xvar [ ]
(∗ : : Subsection : : Closed : : ∗) (∗Support of P l o t t i n g∗)
plotArg [ region_ , v a r l i s t _ : Unevaluated @ xvar [ ] ] :=
Sequence@@Release [ Thread [ Hold [ Prepend ] [ region , Evaluate @ v a r l i s t ] ] ]
SyntaxInformation [ PlotDOA ] = {" ArgumentsPattern " −> {_, _} } ; PlotDOA [ eps_ , region_ ] :=
Block [ {minval , p o i n t s = {} , cp , doa , i n n e r L e v e l S e t } , minval = TimeConstrained [
First@NMinimize [
{MaximalLyapunovFunction , Norm@xvar[]== eps } ,
xvar [ ] , EvaluationMonitor: >(AppendTo[ points , xvar [ ] ] ; ) ] , 1 0 0 ] ;
Print@GraphicsRow [ { ListPlot [ p o i n t s ] , ListPlot [ Function [
Evaluate [ xvar [ ] / . v_[ i_ ]: >ToExpression [ ToString [ v]~~
ToString [ i ] ] ] ,
Evaluate [ MaximalLyapunovFunction / .v_[ i_ ]: >ToExpression [ ToString [ v]~~ToString [ i ] ] ] ] [ # [ [ 1 ] ] , # [ [ 2
] ] ] & / @points ] } ] ;
Print [ minval ] ;
I f [ minval==$Aborted , Return [ ]
I f [ minval <=0,] ; Return [ ] ] ;
(∗
I f [ Length [ FindInstance [ VDotC==0&&VC==minval , xvar [ ] ] ] > 0 , Return [ ]
] ;
∗)
cp = Show[
ContourPlot [ Evaluate@MaximalLyapunovFunction==0, Evaluate@plotArg [ r e g i o n ] , PlotPoints−>40, ContourStyle−>D i r e c t i v e [ Dashed ] ] ,
ContourPlot [ Evaluate@GradientOfMaximalLyapunovFunction==0, Evaluate@plotArg [ r e g i o n ] , PlotPoints−>40,
ContourStyle−>D i r e c t i v e [ Thick ] ] ,
ContourPlot [Norm[ xvar []]== eps , Evaluate@plotArg [ r e g i o n ] , Sequence [ ] (∗ To keep syntax−checker happy ∗) ] ,
AspectRatio−>Automatic , GridLines−>Automatic , GridLinesStyle−>D i r e c t i v e [ Gray, Dotted ] ] ;
doa = RegionPlot [ {
Tooltip [ GradientOfMaximalLyapunovFunction <0, "VDotC" ] , Tooltip [ MaximalLyapunovFunction >=0, "VC" ] ,
Tooltip [Norm[ xvar []] <= eps , "Abs [ x]<=Min@selbord " ]
} , Evaluate@plotArg [ r e g i o n ] , PlotPoints−>40, AspectRatio−>
Automatic , GridLines−>Automatic ,
GridLinesStyle−>D i r e c t i v e [ Gray, Dotted ] , BoundaryStyle−>None ] ;
i n n e r L e v e l S e t = RegionPlot [
Evaluate@GradientOfMaximalLyapunovFunction < 0 &&
Evaluate@MaximalLyapunovFunction < minval && Norm@
L i s t O f V a r i a b l e s <= eps ,
Evaluate@plotArg [ r e g i o n ] , BoundaryStyle−>D i r e c t i v e [ Thick , Dashed ] , GridLines −> Automatic ,
GridLinesStyle −> D i r e c t i v e [ Gray, Dotted ] ] ; Print@Show [
cp ,doa ,
innerLevelSet ,
Graphics [ { PointSize [ Large ] , Point [ { 0 , 0 } ] } ] , GridLines−>
Automatic , GridLinesStyle−>D i r e c t i v e [ Gray, Dotted ] , AspectRatio−>Automatic
] ;
{ eps , minval } ]
(∗ : : Subsection : : ∗)
(∗Functions to f i n d CStar∗)
(∗ checks i f the estimated maximal Lyapunov f u n c t i o n i s not n e g a t i v e at po int p ∗)
maximalLyapunovIsPositiveOrZeroQ [ p_List ] :=
( MaximalLyapunovFunction / . Thread [ xvar []−>p ] ) >= 0
(∗ checks i f g r a d i e n t of the estimated maximal Lyapunov f u n c t i o n i s n e g a t i v e at point p ∗) gradientIsNegativeQ [ p_List ] :=
( GradientOfMaximalLyapunovFunction / . Thread [ xvar []−>p ] ) < 0
(∗ checks i f t h e r e e x i s t p o i n t s i n s i d e the c i r c l e and checkSign f a i l s
∗)
condit [ r_ , checkSign_ ] :=
Quiet [ Reduce [
Exists [ Evaluate@xvar [ ] ,
! checkSign [ xvar [ ] ] && 0 < Norm[ xvar [ ] ] < r ] , Reals
] , {Reduce : : ratnz } ]
Options [ EstimateIntervalOfCStar ] = {" AreaCondition "−>" Condition 1"}
SyntaxInformation [ EstimateIntervalOfCStar ] = {" ArgumentsPattern "−>{_, OptionsPattern [ ] } } ; EstimateIntervalOfCStar [ i n t e r v a l _ L i s t ,
OptionsPattern [ ] ] :=
Module [ {R1value , R2value , midR , midRvalue , temp , fun } , I f [ OptionValue [ " AreaCondition "]==" Condition 1" ,
fun = maximalLyapunovIsPositiveOrZeroQ , fun = gradientIsNegativeQ
] ;
temp = PrintTemporary [ " R1value i s being c a l c u l a t e d . . . " ] ; R1value = ( condit [ i n t e r v a l [ [ 1 ] ] , fun ]=!=False ) ;
NotebookDelete [ temp ] ;
temp = PrintTemporary [ " R2value i s being c a l c u l a t e d . . . " ] ; R2value = ( condit [ i n t e r v a l [ [ 2 ] ] , fun ]=!=False ) ;
NotebookDelete [ temp ] ; Return [
I f [ R1value==R2value , i n t e r v a l ,
( midR = i n t e r v a l [ [ 1 ] ] + ( i n t e r v a l [ [ 2 ] ]−i n t e r v a l [ [ 1 ] ] ) / 2 ;
temp = PrintTemporary [ "midRvalue i s being c a l c u l a t e d . . . " ] ;
midRvalue = ( condit [ midR , fun ]=!=False ) ; NotebookDelete [ temp ] ;
I f [ R1value==midRvalue , {midR , i n t e r v a l [ [ 2 ] ] } ,
{ i n t e r v a l [ [ 1 ] ] , midR}
) ] ] ]
]
(∗ : : Subsection : : Closed : : ∗) (∗Brute−f o r c e check∗)
SyntaxInformation [ ConvergesToTheOriginQ ] = {" ArgumentsPattern " −> {_, _ , _} } ;
ConvergesToTheOriginQ [ solution_ , expectedt_ , epsilon_ ] :=
Module [ {tmax , pmax} ,
(∗ g e t the l e n g t h of the s o l u t i o n ∗) tmax = Min[ s o l u t i o n [ [ All , 2 , 1 , 1 , 2 ] ] ] ;
(∗ i f the s o l u t i o n i s s h o r t e r than expected => considered to be d i v e r g e n t ∗)
I f [ expectedt−tmax>=0.01 expectedt , Return [ False ]
(] ;∗ g e t the endpoint ∗)
pmax = Through [ s o l u t i o n [ [ All , 2 ] ] [ tmax ] ] ; I f [ Norm[ pmax]> e p s i l o n ,
Return [ False ] , Return [ True ] ] ] ;
Options [ ScanDomainOfStability ] = Union [ Options [NDSolve] , {"
ParallelMethod "−>Automatic } ] ;
SyntaxInformation [ ScanDomainOfStability ] = {" ArgumentsPattern "−>{_, _, _, _, _, _, _, _, OptionsPattern [ ] } } ;
ScanDomainOfStability [ system_ , variables_List , t_ , tmax_ , epsilon_ , range_ , resolution_ , solutionIndex_Integer , opts : OptionsPattern [ ] ] :=Module [ {checkOnePoint , grid , g r i d R e s u l t } ,
checkOnePoint [ p0_ ] :=
ConvergesToTheOriginQ [ Quiet@NDSolve [
Evaluate@Join [ system , Thread [ Equal[#[0]&/ @variables , p0 ] ] ] ,
v a r i a b l e s , { t , 0 , tmax } ,
Sequence@@FilterRules [ { opts } , Options [NDSolve ] ] ] [ [ s o l u t i o n I n d e x ] ] ,
tmax , e p s i l o n ] ;
g r i d = Reverse@Transpose [ Outer [ List ,
FindDivisions [ range [ [ 1 ] ] , r e s o l u t i o n [ [ 1 ] ] ] , FindDivisions [ range [ [ 2 ] ] , r e s o l u t i o n [ [ 2 ] ] ] ] ] ;
D i s t r i b u t e D e f i n i t i o n s [ checkOnePoint , ConvergesToTheOriginQ ] ; g r i d R e s u l t = ParallelMap [ checkOnePoint , grid , {2} , Method−>
OptionValue [ " ParallelMethod " ] ] ;
{ gridResult , { g r i d [ [ 1 , {1 , −1}, 1 ] ] , g r i d [[{−1 , 1} , 1 , −1]]} , g r i d }
]
(∗ : : Section : : Closed : : ∗) (∗End of ` Private `∗)
End [ ]
(∗ : : Section : : Closed : : ∗) (∗EndPackage∗)
EndPackage [ ]