# Elegant way to write a system of ODEs with a Matrix

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}];


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?
• Yes you can. But NDSolve can not solve symbolic equations. The coefficient matrix must be numerically. Please correct this for further advice. Commented Jun 1, 2023 at 18:47
• are the ode's all linear? Commented Jun 1, 2023 at 18:48
• @DanielHuber the coefficients l1, l2 etc are specified numerically above the matrix definition in my notebook. @Nasser Yes, all linear. Commented Jun 1, 2023 at 18:49

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.

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]]


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)

• Interesting. How does MatrixExp[] do the computation? For my particular system I am getting an OK solution for t==0 but indeterminate for other times. Commented Jun 1, 2023 at 20:33
• @villaa This works for me: 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]; MatrixExp[ 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] t] . UnitVector[12, 1] ] Commented Jun 1, 2023 at 21:20

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.