11
$\begingroup$

I am trying to numerically solve a time-dependent matrix DE using Mathematica. However, it seems that Mathematica freezes even for two-dimensional matrix equations.

For example, we have:

$$\frac{d}{dt} A(t)=M(t)A(t)+A_0$$

where $M(t)=\begin{pmatrix} 1 & t \\ t & -1 \end{pmatrix}$ and $A_0=A(0) = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$.

For this problem, the simple code below runs forever until Mathematica freezes:

m[t_] = ({{1, t}, {t, -1}}); A0 = ({{1, 0}, {0, 1}});
sA = NDSolveValue[{A'[t] == m[t].A[t] + A0, A[0] == A0}, A, {t, 0, 1}]

I know I can always write it as a $4\times4$ coupled ODEs, but for even slightly higher dimensions the process would be tedious and ugly. Can anyone help?

$\endgroup$
  • 1
    $\begingroup$ The equation in your $\TeX$ expression is not the same as the one you are trying to solve in the code. Trying NDSolveValue[{A'[t] == m[t].A[t], A[0] == A0}, A, {t, 0, 1}] works fine, so which one did you actually need? $\endgroup$ – MarcoB Sep 4 '15 at 21:39
  • $\begingroup$ Thanks for reminding me of the inconsistency. In fact, I want to solve the inhomogeneous time-dependent matrix DE. I've edited my latex equations. $\endgroup$ – K. Qu Sep 5 '15 at 17:23
15
$\begingroup$

The problem is that Dot does not evaluate if one of the arguments is not a list of some sort, and so when the matrix A0 is added, it does something you might not foresee. Consider:

m = {{1, t}, {t, -1}};
A0 = {{1, 0}, {0, 1}};
m.A[t] + A0
(* {
  { 1 + {{1, t}, {t, -1}}.A[t], {{1, t}, {t, -1}}.A[t] }
  , { {{1, t}, {t, -1}}.A[t], 1 + {{1, t}, {t, -1}}.A[t] }
 } *)

which is the same as

(* {
  {1 + m.A[t], m.A[t]}
  , {m.A[t], 1 + m.A[t]}
 } *)

That is clearly not the right expression and is a consequence of the fact that Plus automatically threads over its arguments (as pointed out by Michael E2), whereas Dot remains unevaluated if given symbolic arguments.

Instead, I did this:

m = {{1, t}, {t, -1}};
A0 = {{1, 0}, {0, 1}};
funcs = Array[A[#1, #2][t] &, {2, 2}]
equations = Flatten@Join[
  Thread[D[funcs, t] == m.funcs + A0]
  , Thread[funcs == A0 /. t -> 0]
 ]
sA = NDSolveValue[equations, funcs, {t, 0, 1}]
Plot[sA, {t, 0, 1}]

resulting in

enter image description here

This allows you to generalize to arbitrary dimensions by playing around with the {2, 2} in the definition of funcs. Reconstruct the matrix A[t] by doing

Partition[sA, 2]
$\endgroup$
  • $\begingroup$ Thank you very much! This is exactly what I need. To my understanding, we basically need to treat each matrix element as a function of time and solve the DE. We can use the matrix manipulation to simplify our tedious algebra. This is because NDSolve process does not treat A[t] as a matrix for the matrix operation. Is it correct? $\endgroup$ – K. Qu Sep 5 '15 at 17:33
  • $\begingroup$ @K.Qu. Almost. As pointed out by Marco B in a comment above, without the source term A1, the code works, for reasons basically outlined by Michael E2 in his answer. I tend to get around possible issues like this by solving for the matrix elements rather than the matrix. $\endgroup$ – march Sep 6 '15 at 22:56
4
$\begingroup$

The problem is that Plus will thread symbolic terms, such as m[t].A[t], over explicit arrays, such as A0. One needs all terms to be of the same kind, either all symbolic or all explicit arrays. To code up the symbolic alternative, one can either use DiscreteVariables or code the constant term as a NumericQ-protected constant function.

m[t_] = {{1, t}, {t, -1}};
Clear[A0];
A0 = ({{1, 0}, {0, 1}});
sys = {A'[t] == m[t].A[t] + B[t],
   A[0] == A0, B[0] == A0};
{sA} = NDSolve[sys, A, {t, 0, 1}, DiscreteVariables -> {B}];

Clear[A0];
A0[t_?NumericQ] = {{1, 0}, {0, 1}};
sys = {A'[t] == m[t].A[t] + A0[t],
   A[0] == A0[0]};
{sNA} = NDSolve[sys, A, {t, 0, 1}];

I would have guessed that the DiscreteVariables approach would be easier for NDSolve to analyze and might produce better results, but in this case the results are identical.

sA === sNA
(*  True  *)

With this approach the solution is matrix-valued, which might be an advantage, or not, depending on the applications. It is not an advantage in plotting the components of the matrix:

Plot[A[t] /. sA // Evaluate, {t, 0, 1}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ This is an interesting approach. It forces all the matrices to be explicitly time dependent. I wonder if there is a systematically explanation about the manipulations of matrix or array? I got so confused when a matrix is mixed with functions. $\endgroup$ – K. Qu Sep 5 '15 at 21:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.