# How to speed up a system simulation?

I need to simulate a series of nonlinear (and linear) with thousands of state variables. Here is a simple linear example

{1.00336 x1 + 0.0049723 x2 + 0.00167211 x3, x1, x2}


If I define the following function and then use NestList I get the solution.

model1[{x1_, x2_, x3_}] := {1.003355593149926 x1 +
0.004972295926111435 x2 + 0.0016721109239624015 x3, x1, x2}

NestList[model1, {1, 1, 1}, 2]


The right hand side of model1 is automatically generated by another function. Is there way to do the same with the left hand side? That is, model1 will have the number of arguments depending on the number of x's on the right hand side.

The speed problem arises when the model has thousands of x's. Mathematica does not seem to handle it (abort is the common outcome). Is there a way to handle such models on Mathematica?

Many thanks

• For your first question, you could try defining model1 as model1[x_] := {... x[[1]] + ... x[[2]] + ... x[[3]], ...}, or even model1 = Function[{... Slot[1] + ... Slot[2] + ... Slot[3], ...}] – 2012rcampion Apr 27 '15 at 13:43
• Is the model relatively simple (arithmetic operations only) in all cases? If so, Compile may help you considerably. – Oleksandr R. Apr 27 '15 at 14:01
• @2012rcampion. Many thanks. I think I can modify the generating function as you have suggested. – Ed Mendes Apr 27 '15 at 19:32
• @Oleksandr. Many thanks. Unfortunately it isn't. Could you be so kind to tell me how I could use Compile with model1? Thanks. – Ed Mendes Apr 27 '15 at 19:33
• @2012rcampion. The first option worked ok, although I had to use Quiet to suppress the warning messages. As for the second option, I could not figure out how to do it (it does not give the results as in option 1). Could you so kind to tell me how to use option 2, please? Many thanks – Ed Mendes May 1 '15 at 12:07

Maybe you can try a matrix approach.

1/ The idea is to generate a matrix like this one :

mat = {{a, 1, 0}, {b, 0, 1}, {c, 0, 0}};


then you can see that:

{x1, x2, x3}.{{a, 1, 0}, {b, 0, 1}, {c, 0, 0}}


{a x1 + b x2 + c x3, x1, x2}

gives you the format you want.

2/ Then you can use that matrix directly in NestList (without the need to define a function)

a = 1.003355593149926; b = 0.004972295926111435; c = 0.0016721109239624015;


and

NestList[#.mat &, {1, 1, 1}, 2]


{{1, 1, 1}, {1.01, 1., 1.}, {1.02003, 1.01, 1.}}

returns the same result that your NestList[model1, {1, 1, 1}, 2]

3/ In the case you have hundred a parameters, your matrix will include many identical elements (0's) and it will be probably useful for you to use SparseArray to define the matrix.

Applied to your simple previous example, you have :

matSparse = SparseArray[{{1, 1} -> a, {1, 2} -> 1, {2, 1} -> b,
{2, 3} ->  1, {3, 1} -> c}]


which returns a SparseArray object.

Run Normal[matSparse] if you want to display the full form of the matrix.

You can use this SparseArray object directly in NestList as previously:

NestList[#.matSparse &, {1, 1, 1}, 2]

{{1, 1, 1}, {1.01, 1., 1.}, {1.02003, 1.01, 1.}}


which gives the same result as in 2/.

• It seems ok for a linear system (I shall give it a try). Many thanks. – Ed Mendes Apr 27 '15 at 19:43
• @EdMendes Sorry ;) i focused on your example and "forgot" you need a solution for the nonlinear case. – SquareOne Apr 27 '15 at 19:50