Elegant way to write a system of ODEs with a Matrix

Question

I am using NDSolve[] to solve a first-order system of linear ODEs. The system is a 12-equation system and I have other applications that will be larger. I typically store the coefficients of the ODEs in a matrix so that I can write the system like:

$$ \mathbf{n}^{\prime}(t) = \underline{\underline{\mathbf{m}}} \cdot \mathbf{n}(t)$$

Where $\mathbf{n}(t)$ is a vector of 12 functions and $\underline{\underline{\mathbf{m}}}$ is a 12x12 matrix of coefficients.

To get this into Mathematica I do something like:

m[l1_, l2_, l3_, l4_, l5_, l6_, l7_, l8_, l9_, l10_, l11_, b9a_, 
  b9b_] := -{{l1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-l1, l2, 0, 0, 0,
     0, 0, 0, 0, 0, 0, 0}, {0, -l2, l3, 0, 0, 0, 0, 0, 0, 0, 0, 
    0}, {0, 0, -l3, l4, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, -l4, l5, 0,
     0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, -l5, l6, 0, 0, 0, 0, 0, 0}, {0, 
    0, 0, 0, 0, -l6, l7, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, -l7, l8, 
    0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, -l8, l9, 0, 0, 0}, {0, 0, 0, 0,
     0, 0, 0, 0, -l9*b9b, l10, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 
    0, -l9*b9a, 0, l11, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, -l10, -l11, 
    0}}
    f := {n1[t], n2[t], n3[t], n4[t], n5[t], n6[t], n7[t], n8[t], n9[t], 
  n10[t], n11[t], n12[t]}
    LHS = m[λ1, λ2, λ3, λ4, λ5, λ6, λ7, λ8, λ9, λ10, λ11, b9a, b9b] .f 

Then I follow by creating each equation, doing a command like the one below 12 times.

    Eqns1 = LHS[[1]] == n1'[t] 
    ...

Finally I use the NDSolve[] to get the solution:

sol = NDSolve[{Eqns1, Eqns2, Eqns3, Eqns4, Eqns5, Eqns6, Eqns7, Eqns8,
 Eqns9, Eqns10, Eqns11, Eqns12, n1[0] == 1, n2[0] == 0, 
n3[0] == 0, n4[0] == 0, n5[0] == 0, n6[0] == 0, n7[0] == 0, 
n8[0] == 0, n9[0] == 0, n10[0] == 0, n11[0] == 0, 
n12[0] == 0}, {n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, 
n12}, {t, 0, 1}];

I have a couple of questions about this setup.

1. Can I write the equations more elegantly so that I don't have to write out each equation (12 separate ones in this case)?
2. Is there a way to read the matrix coefficients from a file and make this into a procedure that will work for arbitrary system size?

Answer 1

If you have your system as $x'=Ax+f$ then you can set it up as follows

dep = {x[t], y[t], z[t]};
vars = Head /@ dep ;
mat = {{1, 2, 3}, {0, 0, 3}, {-1, 3, 4}}
fVec = {t, Exp[t], 0}
eqs = Thread[D[dep, t] == mat . dep + fVec ];

ic = {0, 1, 2};
ic = MapThread[#1[0] == #2 &, {vars, ic}]

And now just do

NDSolve[{eqs, ic}, dep, {t, 0, 3}]

So all what you need to do is input the mat and input list of the dependent variables and input the list of initial conditions.

If you want to do the reverse, i.e. you have your equations and want to obtain the matrices, then CoefficientArrays might be useful.

Answer 2

Another alternative, as an explicit matrix-vector system:

m[l1_, l2_, l3_, l4_, l5_, l6_, l7_, l8_, l9_, l10_, l11_, b9a_, 
   b9b_] := -{{l1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-l1, l2, 0, 0, 
     0, 0, 0, 0, 0, 0, 0, 0}, {0, -l2, l3, 0, 0, 0, 0, 0, 0, 0, 0, 
     0}, {0, 0, -l3, l4, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, -l4, l5, 
     0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, -l5, l6, 0, 0, 0, 0, 0, 
     0}, {0, 0, 0, 0, 0, -l6, l7, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
     0, -l7, l8, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, -l8, l9, 0, 0, 
     0}, {0, 0, 0, 0, 0, 0, 0, 0, -l9*b9b, l10, 0, 0}, {0, 0, 0, 0, 0,
      0, 0, 0, -l9*b9a, 0, l11, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 
     0, -l10, -l11, 0}};
Block[{\[Lambda]1, \[Lambda]2, \[Lambda]3, \[Lambda]4, \[Lambda]5, \
\[Lambda]6, \[Lambda]7, \[Lambda]8, \[Lambda]9, \[Lambda]10, \
\[Lambda]11, b9a, b9b},
 {\[Lambda]1, \[Lambda]2, \[Lambda]3, \[Lambda]4, \[Lambda]5, \
\[Lambda]6, \[Lambda]7, \[Lambda]8, \[Lambda]9, \[Lambda]10, \
\[Lambda]11, b9a, b9b} = RandomReal[1, 13];
 sol = NDSolve[{nvec'[t] == 
     m[\[Lambda]1, \[Lambda]2, \[Lambda]3, \[Lambda]4, \[Lambda]5, \
\[Lambda]6, \[Lambda]7, \[Lambda]8, \[Lambda]9, \[Lambda]10, \
\[Lambda]11, b9a, b9b] . nvec[t], 
    nvec[0] == {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}, 
   nvec, {t, 0, 1}]
 ]

(The OP can apologize for all the illegible lambdas.)

The main drawback is that the interpolation returns all 12 vector components. One can extract the (say) 7th component thus

n7Indexed[t_] = Indexed[nvec[t] /. First@sol, 7];

But for a numeric t, it computed all 12 components, then extracts the 7th.

One could reinterpolate if the speed is unsatisfactory. For instance,

n7IFN = Interpolation[
  Transpose@
   Quiet[{nvec["Grid"], nvec["ValuesOnGrid"][[All, 7]], 
      nvec'["ValuesOnGrid"][[All, 7]]} /. First@sol]]

Answer 3

If the matrix is constant numerical, the solution is`

   n[t_]:= Dot[ MatrixExp[t m ], n0 ]

with n0 the list of start values at t=0.

The matrix exponential solves the differential equation

   D[ MatrixExp[t m],t] - m. MatrixExp[t m]  //Simplify  //Chop

It follows, that the constant start vector does not play a role, results can be obtained by Dot[result[t],n0].

While MatrixExp is fine for one value of t, its not suited to be used in a loop like Plot. Use a Module for calculating the exponential as a function of t once.

Module[{mt },mt[t_]=MatrixExp[t m]; Plot

The machinery behind MatrixExp is

     ({diagonal = #1, map = #2} &) @@ Eigensystem[m];

     Total[m - Inverse[map] . DiagonalMatrix[diagonal] .  map, 2]

     mt[t_] = ( Inverse[map] . 
          DiagonalMatrix[Exp[t diagonal]]  . map)

Answer 4

You can solve in vector/matrix form. But you need to inform MA about your data structure with ArrayQ pattern. Try the code below

rhs[f_?ArrayQ, l1_, l2_, l3_, l4_, l5_, l6_, l7_, l8_, l9_, l10_, 
  l11_, b9a_, b9b_] := Module[{m},
   m = -{{l1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-l1, l2, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0}, {0, -l2, l3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 
      0, -l3, l4, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, -l4, l5, 0, 0, 0,
       0, 0, 0, 0}, {0, 0, 0, 0, -l5, l6, 0, 0, 0, 0, 0, 0}, {0, 0, 0,
       0, 0, -l6, l7, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, -l7, l8, 0, 
      0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, -l8, l9, 0, 0, 0}, {0, 0, 0, 0, 
      0, 0, 0, 0, -l9*b9b, l10, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 
      0, -l9*b9a, 0, l11, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, -l10, -l11, 
      0}};
  m.f
  ]
ics = f[0.] == RandomReal[{0., 1.}, 12]
fr = NDSolveValue[{D[f[t], t] == 
    rhs[f[t], 1., 2., 3., 4., 5., 6., 7., 8., 
     9., -0.1, -0.2, -0.3, -0.4], ics}, f, {t, 0., 10.}]
Plot[fr[t][[1]], {t, 0, 10}]

Furthermore, you can introduce some simplifications in the code like pre-construct the matrix.