Solve matrix differential system with mathématica

1. Jul 6, 2011

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

Thanks

2. Jul 6, 2011

Simon_Tyler

Try something like

Code (Text):
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).

Code (Text):

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:
Code (Text):

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

3. Jul 6, 2011

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

4. Jul 6, 2011

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:
Code (Text):

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:
Code (Text):

premier = -Inverse[M].CM;
second = -Inverse[M].KM;
I'll assume that Nul is the zero matrix and ID the identity matrix
Code (Text):

Nul = Array[0 &, {6, 6}];
ID = IdentityMatrix[6];
Code (Text):

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:
Code (Text):
ArrayFlatten[{{premier, second}, {ID, Nul}}]
ArrayFlatten would also be faster since the internal code would be optimized.

5. Jul 6, 2011

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

Code (Text):
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]}]]