Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solve matrix differential system with mathématica

  1. Jul 6, 2011 #1
    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
     
  2. jcsd
  3. Jul 6, 2011 #2
    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}]
     
    kLrxh.gif
     
  4. Jul 6, 2011 #3
    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
     
  5. Jul 6, 2011 #4
    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];
    Then your loop becomes
    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.
     
  6. Jul 6, 2011 #5
    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]}]]
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook