# Need help setting up a numerical simulation and getting output into a table

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

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


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

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

First I solve:

xb[y_]=x/.(Solve[D[f[x,y,sR,nR],x]==0,x]);


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

yC=y/.(Solve[D[g[y,xb[y],sR,nR],y]==0,y]);


Then I want to generate the outcomes

xC=xb[yC]
fC=f[xC,yC,sR,nR]
gC=f[xC,yC,sR,nR]


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 =
{{0.77,0.62,0.502847,0.320859,0.313541,0.191488},
{0.02,0.94,0.266903,0.177778,0.133926,0.0891263}}


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.

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} =
Solve[
{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; 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}}


### Update

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; 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}}

• 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 .... – Bjoj Jan 11 '17 at 12:38
• ...Also, to understand exactly how this works, do you mind explaining what is the function of [[1, All, 2]] in your code? – Bjoj Jan 11 '17 at 13:16
• @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. – m_goldberg Jan 11 '17 at 22:08
• Thank you! This was very informative – Bjoj Jan 13 '17 at 13:31
• 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? – Bjoj Jan 16 '17 at 8:50