# DSolve for large nearly sparse matrix

The following code creates a system of first order linear homogeneous differential equations and attempts to solve them. The time needed scales as the 4th power of the matrix size. Is there an approach to solving this problem that would be more efficient? Thanks!

itest = 10;
h = Table[i, {i, itest}];

t0 = AbsoluteTime[];
bmtest = Join[
Table[RandomReal[{0, .25}]
, {i, itest - 1}, {j, itest}]
, {Join[
Table[0, {i, itest - 1}]
, {-1}]
}];
For[i = 1, i <= itest - 1, i++,
bmtest[[i, -1]] *= 10;
bmtest[[i, i]] *= -100;
av = Complement[h, {i, itest}];
rs = Join[RandomSample[av, 2], {i}];
For[j = 1, j <= itest - 1, j++,
If[! MemberQ[rs, j],
bmtest[[i, j]] = 0]
]
];

t1 = AbsoluteTime[];
eqn = Table[{D[Subscript[n, i][t], t]}, {i, 1, itest}] ==
bmtest.Table[{Subscript[n, i][t]}, {i, 1, itest}];
ic = Join[
Table[Subscript[n, i][0] == 0, {i, 1,
itest - 1}], {Subscript[n, itest]'[0] == -1}];
vars = Table[Subscript[n, i][t], {i, 1, itest}];
dsv = DSolve[{eqn, ic}, vars, t];
t2 = AbsoluteTime[];

Plot[Evaluate[Table[dsv[[1, i, 2]], {i, 1, itest}]], {t, 0, 6}
, PlotStyle -> Table[Hue[i/itest], {i, 1, itest}]
, PlotRange -> All]

• "system of first order linear homogeneous differential equations" - have you looked at MatrixExp[] already? Plot[Evaluate[MatrixExp[bmtest t, UnitVector[itest, itest]]], {t, 0, 6}, PlotStyle -> Table[Hue[i/itest], {i, 1, itest}], PlotRange -> All] Dec 11 '21 at 9:45
• @J.M. I have not looked into MatrixExp[]. I will read about it and get back to you. Thank you for the suggestion. Dec 11 '21 at 18:21
• @J.M. A little bit of testing shows that MatrixExp[] scales as the 3rd power of the matrix dimensions, vs the 4th power for DSolve[]. This is exponentially better (literally!). 6.5 sec to solve 100 differential equations, 52.5 sec to solve 200. If the scaling is consistent, I should be able to solve 500 differential equations (the goal) in 14 minutes. Dec 12 '21 at 20:59
• That sounds good. The only thing you need to do now is to actually use SparseArray[] to represent your matrix, instead of as a list. Dec 13 '21 at 1:08
• @J.M. I turned my system into a SparseMatrix. I'm finding that the time needed to solve the system of equations as a SparseMatrix or as a List is the same. They both scale as the 3rd power of the matrix dimensions, as long as I'm using MatrixExp[]. Perhaps this is the best that can be done (?). Dec 16 '21 at 17:58