Solve matrix differential system with mathématica

[Nesrine]

hello ,

I need to solve a differential system with mathematica that has the form :

B'=A*B

knowing that B and A are (12,12) Matrix and with initial conditions B(0)=I

Can you help me please

Thanks

 

[Simon_Tyler]

Try something like

n = 12;
B[t_] = Array[Subscript[b, ##][t] &, {n, n}];
Bt = Array[Subscript[b, ##] &, {n, n}];
id = IdentityMatrix[n];
(*A=2.id-RotateRight[id,1]-RotateRight[id,-1];*)
A = SparseArray[{Band[{1, 1}] -> 2., Band[{1, 2}] -> -1, Band[{2, 1}] -> -1}, {n, n}];

The key command for getting a matrix equation into something that DSolve and friends like is the command LogicalExpand. (You can also use combinations of Flatten and Thread, but that's not as neat).

soln = DSolve[LogicalExpand[B'[t] == A.B[t] && B[0] == id], Flatten@Bt, t];

Where we could have used B[t] instead of Bt, but then the resultant replacement rules (soln) would be less useful. For the purely numerical example above, NDSolve would be quicker to get a solution for any particular range of t.

We can then plot the result using ArrayPlot:

Animate[ArrayPlot[Evaluate[B[t] /. soln[[1]]]], {t, 0, 5 n}]

 

[Nesrine]

Thank you very much for your help...but the problem is that since I'm not used to Mathematica I don't understand so functions that you have used.

For example I supposed in this part you defined the Matrix A
(*A=2.id-RotateRight[id,1]-RotateRight[id,-1];*)
A = SparseArray[{Band[{1, 1}] -> 2., Band[{1, 2}] -> -1, Band[{2, 1}] -> -1}, {n, n}];

But for me the Matrix A has been calcultaed before.
Can you show me please how to implemented in my programm solution that you proposed to me .

Here is my programm for calcultaing the matrix A :

A = Table[0, {12}, {12}];
premier = -Inverse[M] . CM;
second = -Inverse[M]. KM;
For[k = 1, k <= 3, k++,
For[l = 1, l <= 3, l++,
A[[k]][[l]] = premier[[k]][[l]];
A[[k]][[l + taille]] = second[[k]][[l]];
A[[k + taille]][[l]] =
ID[[k]][[l]]; (*Calcul de la matrice d'état*)


A[[k + taille]][[l + taille]] = Nul[[k]][[l]]
]
];

Where M, CMand KM are periodic matrix (6,6)

Thank you very much I really appreciate

 

[Simon_Tyler]

I only defined A like that since I needed some coefficient matrix. You can see what I (almost) arbitrarily chose by running MatrixForm[A].
By the way, the stuff in (* ... *) is a comment - an alternative choice of A that I should have deleted before posting.

So to make the code I posted in your program work, you just need your own A.
There are quite a few things not defined in your code, which makes it hard to tell what you want.
It kind of looks like you're constructing the block matrix
{{premier, second}, {ID, Nul}}
but then your iterators should go up to 6...
Anyway, I'll assume that's what you're doing.
If it's not, can you describe in words how you construct A and what the result should be?

Let's first construct some test cyclic matrices:

M  = NestList[RotateRight, RandomReal[{0, 1}, 6], 5];
CM = NestList[RotateRight, RandomReal[{0, 1}, 6], 5];
KM = NestList[RotateRight, RandomReal[{0, 1}, 6], 5];

and your premier and second matrices:

premier = -Inverse[M].CM;
second = -Inverse[M].KM;

I'll assume that Nul is the zero matrix and ID the identity matrix

Nul = Array[0 &, {6, 6}];
ID = IdentityMatrix[6];

Then your loop becomes

A = Array[0 &, {12, 12}];
taille = 6;
For[k = 1, k <= taille, k++,
  For[l = 1, l <= taille, l++,
   A[[k]][[l]] = premier[[k]][[l]];
   A[[k]][[l + taille]] = second[[k]][[l]];
   A[[k + taille]][[l]] = ID[[k]][[l]];
   A[[k + taille]][[l + taille]] = Nul[[k]][[l]]
   ]];

Which is much more easily created using ArrayFlatten:

ArrayFlatten[{{premier, second}, {ID, Nul}}]

ArrayFlatten would also be faster since the internal code would be optimized.

 

[Simon_Tyler]

Oops... I confused periodic matrix (M^(k+1) = M) with cyclic matrix (M_{i,j}=M_{i+1,j+1}) - not that it makes any difference to the above example.

Also, if Nul really is the zero matrix then you can replace Nul[[k]][[l]] with 0 in the loop and Nul with 0 in the ArrayFlatten.

Also, when accessing elements of arrays, it is about twice as fast to use M[[k,l]] than M[[k]][[l]]. To check this, run

M = RandomReal[{0,1},{6,6}]
Timing[Do[M[[k]][[l]],{k,RandomInteger[{1,6},1000]},{l,RandomInteger[{1,6},1000]}]]
Timing[Do[M[[k,l]],{k,RandomInteger[{1,6},1000]},{l,RandomInteger[{1,6},1000]}]]