# How do I represent a system dynamics feedback loop?

In System Dynamics, if I want represent the relationship between Speed and Distances, I create a Flow (Speed) and a Stock (Distance) as you can see in this Insight Maker Sample . Heres an Image of how it looks like:

Now, keeping a constant speed of 20 results in this nice Graph:

I know how to represent that in Mathematica, it is easy:

speed = 20
distance = Integrate[speed, time]
Plot[distance, {time, 0, 20}]


That will get me this nice graph:

Now I want to do something more complex, I want to add a feedback loop, I have already done it in InsightMaker too, it looks like this:

• The flow Reproduction is configured as: 20*[Reproduction Rate]
• The stock E. Coli Bacteria is configured as: Initial Value: 1
• The variable Reproductino Rate is configured as: [ E. Coli Bacteria]

So I get exponential growth :

My question is: how do I represent this kind of system dynamics model with a feedback loop in Mathematica?

UPDATE: Here is a table with the values generated by InsightMaker

• Hour 0 E. Coli Bacteria: 1
• Hour 1 E. Coli Bacteria: 21
• Hour 2 E. Coli Bacteria: 441
• Hour 3 E. Coli Bacteria: 9,261
• Hour 4 E. Coli Bacteria: 194,481

Would you be so kind as to please configure your answers to match this results? (I am learning as a hobby system dynamics & mathematica and my goal with this question is to understand how to translate a system dynamics diagram to mathematica, so I want to be sure I correctly understad how each part of a system dynamics diagram matches an equivalent mathematica code)

-

Of course the true exponential function is the solution of a differential equation in which the time steps are taken to be infinitesimal. However, it seems like you're looking for an implementation that's as literal as possible and uses discrete time steps:

coli[rate_, nStart_, n_] :=
NestList[Function[{coli}, Round[coli + rate coli]], nStart, n]

ListLinePlot[coli[1, 1, 10]]


This is just a rough example, the main ingredient is NestList which repeatedly applies a function to a starting configuration. This is a way of realizing the "loop" described by the graphic you showed. Above, nStart is the initial number of bacteria, and n the number of time steps. Since the step size is unity, the rate is measured in a dimensionless way relative to the time step.

Edit

It is a trivial matter to adjust the above to the given parameters that have been posted in the update to the question:

coli[20, 1, 4]

(* ==> {1, 21, 441, 9261, 194481} *)


Edit 2

In response to the question about the time variable, one could add some additional flexibility to the function that is being nested here, by keeping track of a time variable explicitly that serves as a clock. In this case it's simply incremented by 1 at every step. Then the starting configuration of the loop has to be a tuple {startTime, startNumberoFColi} = {0, 1}:

coli[rate_, nStart_, n_] :=
NestList[
Function[
{timeColi},
{
timeColi[[1]] + 1, Round[timeColi[[2]] + rate timeColi[[2]]]
}
],
{0, nStart}, n]

coli[20, 1, 4]

(* ==> {{0, 1}, {1, 21}, {2, 441}, {3, 9261}, {4, 194481}} *)

ListLinePlot[%]


The Function now takes the tuple of time and number as its argument timeColi, and I have to extract the time from it by saying timeColi[[1]], while the number is the part timeColi[[2]]. They are then re-calculated individually and returned as a tuple, ready for the next iteration.

-
I like it... But the number of bacteria in your graph does not match mine... –  Luxspes Apr 6 '13 at 13:03
You are right when you say that I am looking for an implementation that's as literal as possible, the better the solution matches the diagram parts, the better it will help me understand the math behind the diagram –  Luxspes Apr 6 '13 at 17:34
All you need to do is choose the rate and starting value to be the same as in your question, and use the function I posted. I've updated the answer. –  Jens Apr 6 '13 at 17:41
I tried with coli[20, 1, 4], and the numbers match! {1, 21, 441, 9261, 194481} but there is a mismatch with the time, 194481 is for hour 5, but in InsightMaker is for hour 4... –  Luxspes Apr 6 '13 at 17:43
Maybe that's because I'm in a different time zone? No, the counter variable n in NestList simply starts at 1 while your clock starts at 0. It's just a numbering convention. –  Jens Apr 6 '13 at 17:52

This might not be the canonical way, but the way I would do it is to represent the system as a pair of state-updating equations, translate that into matrix form and use NestList to show the time path. This is of course assuming that you are happy to work in discrete time.

Consider, for example, where the reproduction rate at time $t$ is a positive function of the rate at time $t-1$ and a negative function of the total population at time $t-1$.

$\Delta p_t = \beta_1 \Delta p_{t-1} - \beta_2 p_{t-1}$

Where $\beta_1$ and $\beta_2$ are positive numbers.

The population itself is a simple stock-updating equation.

$p_t = \Delta p_{t-1} + p_{t-1}$

Translating this to Mathematica syntax implies a matrix of coefficients

{{beta1, -beta2},{1,1}}


Say, for example:

coefs = {{0.9, -0.02}, {1, 1}};


Then you can make this into a simple NestList like this:

NestList[coefs.#1 &, {0.01, 1.}, 20]


However the result for this particular configuration goes negative, suggesting that some additional logic is needed:

result = NestList[With[{rawnext = coefs.#1}, {First@rawnext, Max[Last@rawnext, 0]}] &,
{0.01, 1.}, 20]


You can plot the population by taking the second column of the results.

ListLinePlot[result[[All, 2]] ]


But what's really nice about Mathematica is that you can explore the parameter space very easily, using Manipulate.

Manipulate[
ListLinePlot[
NestList[With[{rawnext = {{beta1, beta2}, {1, 1}}.#1}, {First@
rawnext, Max[Last@rawnext, 0]}] &, {delta0, pop0}, 20][[All, 2]] ],
{{beta1, 0.99}, 0.01, 2}, {{beta2, 0.01}, -1, 2}, {delta0, 0, 1}, {pop0, 0.5, 10}]


Another way do simulate systems is to use things like RSolve. But for simple cases like this, a matrix representation is easy to implement.

-
Thanks!, but one thing bothers me, after 4 hours the number of bacteria in your graph does not match mine... –  Luxspes Apr 6 '13 at 13:01
I am not exactly sure of how to translate "beta1", "beta2" , "delta0" and "pop0" to the diagram parts... could you please change the name of those parameters to better match the diagram parts? –  Luxspes Apr 6 '13 at 17:33