4
$\begingroup$

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

4 Answers 4

7
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

And now just do

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

Mathematica graphics

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.

$\endgroup$
5
$\begingroup$

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

Mathematica graphics

(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]]
$\endgroup$
1
$\begingroup$

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)
$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – villaa
    Jun 1, 2023 at 20:33
  • $\begingroup$ @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] ] $\endgroup$
    – Michael E2
    Jun 1, 2023 at 21:20
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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