# How to initiate a large system of differential equations?

I want to numerically solve a system of n differential equations of the form

x''[t] == M.x[t] + F[t], x' == Table[0,{i,1,n}], x == Table[RandomReal[{-0.1,0.1}],{i,1,n}]


where M is some n times n matrix. If f[t]=0 everything works just fine, but if I put in for example (n = 3)

F[t] = {f[t],0,0}


nothing works and I get memory overflow. Now if I write down each equation separately it works just fine

x1''[t] == -2*x1[t] + x2[t] + x3[t] + f[t]
x2''[t] == -2*x2[t] + x1[t] + x3[t]
x3''[t] == -2*x3[t] + x1[t] + x2[t]


The problem is that I want to use quite large n and different matrices M. Therefore it is inconvenient to write down all equations by hand. Is there a more compact way of writing down those equations with the inhomogeneity at one site?

• Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Nov 23 '16 at 14:59
• I think to get useful answers you would need to give more details. When I try your example with M as read from your 3-equations-example, with n=3 and f[t_]:=0.1*Sin[t] NDSolve solves with no errors or problems for 0<=t<=1. That is admittedly a very simple case but it shows that basically what you try should work and that only for your specific settings there seems to be a problem -- without those we can't help. We would at least nee to know: definition of M, definition of f, value of n, time range for which you try to solve for. – Albert Retey Nov 23 '16 at 17:48
• There are a few examples oft this kind in the Finite Element Programming tutorial. – user21 Nov 24 '16 at 7:19

I don't understand what you specifically are looking for, but is it something like this? Careful, I changed the system to a first-order system (I assumed that you know that you can always do that).

n = 5;
M = RandomReal[{-1, 1}, {n, n}];
x0 = RandomReal[{-1, 1}, n];
f = RandomSample[Table[t^p, {p, 0, 10}], n];
xf = Table[ToExpression["x" <> ToString[i] <> "[t]"], {i, n}];
eqs = D[xf, {t, 1}] - (M.xf + f);
eqs = Table[eqs[[i]] == 0, {i, Length@eqs}];
cond = Table[(xf[[i]] /. t -> 0) == x0[[i]], {i, Length@x0}];
xsol = NDSolveValue[Flatten[{eqs, cond}], Evaluate@xf, {t, 0, 1}];
Plot[Evaluate[xsol /. t -> tp], {tp, 0, 1}, PlotLegends -> Automatic] • You could use Thread[eqs == 0] to avoid one or two Table. – anderstood Nov 23 '16 at 17:10
• Cool, thank you very much! – Mauricio Fernández Nov 23 '16 at 17:42
• I think it would be worth explaining that what you do is to construct the system of equation programmatically -- something that as a symbolic algorithm is very powerful and somewhat specific to Mathematica. It might be a valid workaround if there really is a problem with what R. Milus tries to do, but I think what he tries should also be possible with his original approach (solving one eq. for a vector-valued function)... – Albert Retey Nov 23 '16 at 17:55

Here is an example which shows that in principle what you try to do should work:

n = 20;
M = SparseArray[
{{i_, i_} -> -2., {i_, j_} /; Abs[i-j] == 1 -> 1.}, {n,n}, 0.
];
x0 = RandomReal[{-1, 1}, n];
F[t_] = RandomChoice[{0.1*Sin[t], 0}, n];

xsol = NDSolveValue[{
x''[t] == M.x[t] + F[t],
x' == ConstantArray[0, n],
x == x0
}, x, {t, 0, 1}
];
Plot[xsol[t], {t, 0, 1}]


Of course depending on values/definitions of n, M and f there still could be problems but I would suspect that those then would depend on other things than whether you formulate the problem as one equation for one vector-valued function or as many equations for many scalar-valued functions. For me the system doesn't look very well behaved with these settings as n increases...

Mauricio Lobos has shown in his answer that it is also a valid and good approach in Mathematica to construct an array of equations programmatically which is possible due to the symbolic nature of the language. It probably is worth noting that since several versions of Mathematica it is not necessary that each dependent variable is a symbol anymore, so when creating such equations programmatically it is now much easier (and does away with several problems that generating symbols on the fly has) to not create a single symbol with ToExpression for every equation but instead use something like:

depvars = Table[x[k], {k, n}];
equations = Join[
Thread[D[Through[depvars[t]], {t, 2}] == (M.Through[depvars[t]] + F[t])],
Thread[Derivative[#] & /@ depvars == 0],

On the other hand there could be large differences in the runtimes and memory-efficiency of the two approaches. But actually I would expect that for large n the vector-valued function approach would have the potential to be much more efficient, especially concerning memory-usage. But that depends on implementation details that I don't have insight to and I have not done any experiments about that, so you might want to try out which approach will work better for you.