I would like to set up a simple numerical simulation of a model - but I don't know how to do it in a simple and efficient manner. Here's what I want to do (using a simplified example):

I want to first draw random values for the different parameters in my model, here I assume a uniform distribution between zero and one:

sR = RandomReal[{0, 1}];
nR = RandomReal[{0, 1}];

Next I want to define some objective functions ...


and then solve a some equations, using the functions I just defined:

(EDIT: Originally my example involved solving the equations simultaneously)

First I solve:


Then I want to solve the next equation, substituting in the solution I got from the first:


Then I want to generate the outcomes


I want to run this operation multiple times and put all the data points in a matrix, with each row showing {sR,nR,xC,yC,fC,gC}: As an example, suppose I want to repeat the operation two times, the first time I draw the values $sR=0.77$ and $nR=0.62$, and the second time I draw $sR=0.02$ and $nR=0.94$. In this case I want the final data to look like this

m = 

I know one way of getting the output that I want - but it involves defining the functions and deriving the solutions on the more general form first, before generating a random dataset for the parameter values, and using the MapThread function to generate the final dataset. This method can be very inefficient (and eats a lot of memory) when the functions and equations are a lot more complex - and when I want to generate a large dataset of say 1k-200k rows. For the model that I want to run, Mathematica did not manage to finish before my computer was out of memory.

Maybe an alternative method is to use the Table function(?) - but I don't know how to do it in my example.


1 Answer 1


It seems simple to me. But we need make changes in your definitions.

f[x_, y_, sR_, nR_] := x (1 - 2 nR*x + sR*y)
g[y_, x_, sR_, nR_] := y (1 - 3 nR*y + .5 sR*x)

simStep :=
  Module[{sR, nR, xC, yC},
    {sR, nR} = RandomReal[1, 2];
    {xC, yC} =
        {D[f[x, y, sR, nR], x] == 0, 
         D[g[y, x, sR, nR], y] == 0}, 
        {x, y}][[1, All, 2]];
    {sR, nR, xC, yC, f[sR, nR, xC, yC], g[sR, nR, xC, yC]}]

Now we can use Table to run your simulation.

SeedRandom[1]; Table[simStep, 6]
{{0.817389, 0.11142, -41.1394, -23.6545, 28.679, 46.3565}, 
 {0.789526, 0.187803, 3.58326, 2.14279, -1.35059, -2.95196}, 
 {0.241361, 0.0657388, 8.52379, 5.14322, -0.222632, -0.589874}, 
 {0.542247, 0.231155, 1.69917, 1.05318, 0.135889, -0.280269}, 
 {0.396006, 0.700474, 0.393148, 0.256456, 0.424627, 0.329881}, 
 {0.211826, 0.748657, 0.350263, 0.23088, 0.246653, 0.20852}}


The simulation model as reformulated by the OP can be implemented using the same approach as used for the original model, but we must be more careful about creating the local scopes for the variables. In this case, some variables need to be manipulated symbolically; they are not just holders of numeric values, so we need to localize them with Block.

altSimStep :=
  Block[{x, y, xb},
    Module[{sR, nR, xC, yC},
      {sR, nR} = RandomReal[1, 2];
      xb[y_] = Solve[D[f[x, y, sR, nR], x] == 0, x][[1, 1, 2]];
      yC = Solve[D[g[y, xb[y], sR, nR], y] == 0, y][[1, 1, 2]];
      xC = xb[yC];
      {sR, nR, xC, yC, f[sR, nR, xC, yC], g[sR, nR, xC, yC]}]]

SeedRandom[1]; Table[altSimStep, 6]
{{0.817389, 0.11142, -1.98915, -2.30798, 3.72028, 5.35288}, 
 {0.789526, 0.187803, 6.72905, 5.13593, -4.6157, -8.31606}, 
 {0.241361, 0.0657388, 11.5484, 8.43843, -0.558567, -1.14177}, 
 {0.542247, 0.231155, 1.79104, 1.20984, 0.0552789, -0.412699}, 
 {0.396006, 0.700474, 0.393392, 0.258187, 0.424152, 0.329101}, 
 {0.211826, 0.748657, 0.35029, 0.231266, 0.246623, 0.20847}}
  • $\begingroup$ Thank you! This works - and I think I can make it work also in my case. But I forgot something in my original question, if you don't mind looking at my edit? What if I want to solve the equations in steps, using the output from the first one when solving the second one? I tried to modify your code, but it didn't work .... $\endgroup$
    – bjoj
    Jan 11, 2017 at 12:38
  • $\begingroup$ ...Also, to understand exactly how this works, do you mind explaining what is the function of [[1, All, 2]] in your code? $\endgroup$
    – bjoj
    Jan 11, 2017 at 13:16
  • $\begingroup$ @Bjoj. 1) I have updated my answer to address the reformulated simulation model. 2) [[1 All, 2]] is a efficient way to extract the values from the list of rules returned by Solve. Look up Part in the documentation. $\endgroup$
    – m_goldberg
    Jan 11, 2017 at 22:08
  • $\begingroup$ Thank you! This was very informative $\endgroup$
    – bjoj
    Jan 13, 2017 at 13:31
  • $\begingroup$ Do you mind if I ask you one final question? I ran a few tests using the solution you provided, and I found that it doesn't seem to matter at all whether I localize the variables in Block or not. Specifically, I get the exact same output whatever I use of the following three options at the start: 1) Module[{sR,nR,xb,yC},, 2) Module[sR,nR,xb,yC,xC] or 3) Block[{x,y,xb},Module[{sR,nR,xC,yC}, (the solution you provided). Can you explain when and why it matters, and perhaps why it doesn't seem to matter in my case? $\endgroup$
    – bjoj
    Jan 16, 2017 at 8:50

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.