How to speed up a system simulation?

Question

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

Answer 1

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)

In your case :

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/.