- #1

Nesrine

- 9

- 0

I need to create a program on mathématica 8 to study the stability of my system using Floquet transition matrix .

to compute the Floquet transition matrix I made this program based on a fourth order Runge Kutta integration :

X1[t_] = {x1[t], x2[t], x3[t], x4[t], x5[t], x6[t] , x7[t], x8[t],

x9[t], x10[t], x11[t], x12[t]};

system1 = MapThread[#1 == #2 &, {X1'[t], A.X1[t]}];

X2[t_] = {x1[t], x2[t], x3[t], x4[t], x5[t], x6[t] , x7[t], x8[t],

x9[t], x10[t], x11[t], x12[t]};

system2 = MapThread[#1 == #2 &, {X2'[t], A.X2[t]}];

X3[t_] = {x1[t], x2[t], x3[t], x4[t], x5[t], x6[t] , x7[t], x8[t],

x9[t], x10[t], x11[t], x12[t]};

system3 = MapThread[#1 == #2 &, {X3'[t], A.X3[t]}];

X4[t_] = {x1[t], x2[t], x3[t], x4[t], x5[t], x6[t] , x7[t], x8[t],

x9[t], x10[t], x11[t], x12[t]};

system4 = MapThread[#1 == #2 &, {X4'[t], A.X4[t]}];

CRK4[]["Step"[rhs_, t_, h_, y_, yp_]] :=

Module[{k0, k1, k2, k3}, k0 = h yp;

k1 = h rhs[t + h/2, y + k0/2];

k2 = h rhs[t + h/2, y + k1/2];

k3 = h rhs[t + h, y + k2];

{h, (k0 + 2 k1 + 2 k2 + k3)/6}]

CRK4[___]["DifferenceOrder"] := 4

dstep1 = NDSolve[{system1, x1[0] == 1,

x2[0] == x3[0] == x4[0] == x5[0] == x6[0] == x7[0] == x8[0] ==

x9[0] == x10[0] == x11[0] == x12[0] == 0}, {x1, x2, x3, x4, x5,

x6, x7, x8, x9, x10, x11, x12}, {t, 0, 2 \[Pi]},

Method -> {"DoubleStep", Method -> CRK4}];

dstep2 = NDSolve[{system2, x2[0] == 1,

x1[0] == x3[0] == x4[0] == x5[0] == x6[0] == x7[0] == x8[0] ==

x9[0] == x10[0] == x11[0] == x12[0] == 0}, {x1, x2, x3, x4, x5,

x6, x7, x8, x9, x10, x11, x12}, {t, 0, 2 \[Pi]},

Method -> {"DoubleStep", Method -> CRK4}];

dstep3 = NDSolve[{system3, x3[0] == 1,

x1[0] == x2[0] == x4[0] == x5[0] == x6[0] == x7[0] == x8[0] ==

x9[0] == x10[0] == x11[0] == x12[0] == 0}, {x1, x2, x3, x4, x5,

x6, x7, x8, x9, x10, x11, x12}, {t, 0, 2 \[Pi]},

Method -> {"DoubleStep", Method -> CRK4}];

dstep4 = NDSolve[{system4, x4[0] == 1,

x1[0] == x2[0] == x3[0] == x5[0] == x6[0] == x7[0] == x8[0] ==

x9[0] == x10[0] == x11[0] == x12[0] == 0}, {x1, x2, x3, x4, x5,

x6, x7, x8, x9, x10, x11, x12}, {t, 0, 2 \[Pi]},

Method -> {"DoubleStep", Method -> CRK4}];

and after that I get my transition matrix by using the solutions as the column.

The problem is that I don't get the expected result , so I'm wondering if my understanding of the method is correct or not.

If someone know how to compute Floquet transition matrix in mathematica 8 , I'll be very greatfull

thxx