On and off over the last few months I've been meaning to learn about computing a center manifold approximation and normal form of a dynamic system of three differential equations.

My main reading: Stephen Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Second Edition, Springer, 2003.

I want to apply the technique to a system I derived from an optimal control problem. As a first step, I decided to reproduce the steps for the following system, for which the solution has been published.

ydot := diff(y(t),t) = -y(t)*z(t) -s*y(t)*q(t);
qdot := diff(q(t),t) = (1+q(t))*(y(t)-s*q(t));
zdot := diff(z(t),t) = s*q(t)+z(t)*(y(t)-s*q(t));
sys:=[ydot,qdot,zdot];

The article I have at hand (reference upon request) applies the following transformation.


[ y ]            [ s   0    0 ]       [ u ]
[   ]            [            ]       [   ]
[ q ]     =      [ 1    0   1 ]       [ v ]
[   ]            [            ]       [   ]
[ z ]            [ 0   1   -1 ]       [ w ]

The system may then be transformed as:

 

[ u' ]         [ 0   0    0 ]      [ u ]

[    ]         [            ]      [   ]

[ v' ]    =    [ s    0   0 ]      [ v ]   +  Higher Order Terms

[    ]         [            ]      [   ]

[ w' ]         [ 0    0  -s ]      [ w ]

I would like to obtain that transformation with Maple. In the following I derived the Jordan normal form, but it isn't the same. I just don't know how to get the above transformation. I'll be working on that over the next few days.


F := subs(y(t)=y,q(t)=q,z(t)=z,map(rhs,sys));
J := VectorCalculus[Jacobian](F,[y,q,z]);
J0 := eval(J,[y=0,q=0,z=0]);
E0,V0 := LinearAlgebra:-Eigenvectors(J0);
Jjordan,Q := LinearAlgebra:-JordanForm(J0,output=['J','Q'] );
'Q^(-1)' = Q^(-1);
''Jjordan'' = Q^(-1) . J0 . Q: convert(%,rational);  #check: returns Jjordan
''J0'' = Q . Jjordan . Q^(-1): convert(%,rational); #check: returns J0

      F := [-y z - s y q, (1 + q) (y - s q), s q + z (y - s q)]


                               [0    0     0]
                               [            ]
                         J0 := [1    -s    0]
                               [            ]
                               [0    s     0]


                              [0 ]  [0    0     0]
                              [  ]  [            ]
                    E0, V0 := [0 ], [0    0    -1]
                              [  ]  [            ]
                              [-s]  [1    0     1]


                       [0    1    0 ]  [0      1        0  ]
                       [            ]  [                   ]
         Jjordan, Q := [0    0    0 ], [0     1/s     - 1/s]
                       [            ]  [                   ]
                       [0    0    -s]  [1    - 1/s     1/s ]


 


Please Wait...