Stocks, flows, and integration

Manual contents | SDSir tutorial | SDMain reference | Stock and Store reference

This chapter maps classic system-dynamics concepts onto their NuWebKit code constructs, explains exactly what happens in one tick of the clock, and shows how the same machinery serves continuous (Euler/RK4) and discrete/stochastic models. Sources: sdmain.js, stock.js, clock.js, integrator.js.

1. From diagram to code

Two consequences of this mapping are worth internalizing:

• A flow is a function, not a number. The runtime calls it whenever the integrator needs the rate — once per dt under Euler, four times per dt under RK4. It must therefore be a pure computation over stock levels and parameters.
• A biflow or shared flow appears in every stock it touches. SD tools draw one pipe between two stocks; in code the same term occurs twice, negative in the source and positive in the sink. Conservation checks (do the stocks sum to a constant?) catch sign and omission errors.

2. Stocks, auxiliaries, and parameters

A Stock holds a two-phase value: the committed current level that get() returns, and a pending next contribution the integrator fills in. During one tick every flow reads the same committed levels; new levels are committed together afterward. That is what makes the update synchronous — the order in which stocks were created never changes the answer. The machinery is documented in the Stock and Store reference.

Auxiliaries need no special class. In the Lotka-Volterra model below, predation() is an auxiliary used by two flows; defining it once keeps the diagram's structure visible in code. Auxiliaries recompute from committed state whenever called, which is exactly the SD semantics of a converter.

Parameters arrive from the Project's inputs, keyed by code. Read them through getters so every flow evaluation sees the latest slider position; copying a parameter into a field at init() time would freeze it.

3. The simulation cycle

Under Euler (visitMode: "SD"), each clock tick performs, in order:

1. step() — your per-tick hook. The base class calls coordinate() from here in Euler mode; both exist for discrete logic (see section 7). For purely continuous models they do nothing.
2. Flow evaluation. Every stock computes its pending change: next = dt * flow(). All flows read the same committed levels from get(); nothing has moved yet.
3. Commitment (shift). Every touched stock commits: current = current + next. This is the forward-Euler update x(t + dt) = x(t) + dt * f(x(t)), applied to all stocks simultaneously.
4. Time advances by dt.
5. update() — your output hook: read the new levels and show them.
6. The runtime sends buffered display data to the charts and checks the stop conditions (maxTime, stopWhenTrue).

Under RK4 the same skeleton holds, but phase 2–3 becomes the four-stage Runge-Kutta sequence described in section 5; coordinate() is not called in RK4 mode, and update() still runs once per completed dt.

The practical rules this ordering implies:

• Inside a flow, get() always returns the level at the start of the current step (or RK4 stage) — never a half-updated value.
• update() is the first place the new levels are visible; it is the right place for show, logging, and derived observations.
• Assigning state or producing output inside a flow breaks the contract: the integrator may call it four times per dt, or your side effect lands between stages.

4. Time step and accuracy

dt scales every Euler step and every RK4 stage. Project dt sets it; the hosting page can change it live (the runtime re-installs it into every stock), and maxTime sets the stop time — the run halts once the clock reaches it, within half a dt.

Choosing dt is the classic SD trade-off: too large and the discretization distorts the dynamics (oscillations flatten, peaks shift, stocks can overshoot negative); too small and runs crawl. The standard check applies unchanged: halve dt and rerun — if trajectories change visibly, dt was too large. The SDSir tutorial experiments demonstrate this concretely.

5. Euler and RK4 integration

visitMode selects the integrator for the whole model:

Under RK4 the runtime evaluates every flow at the step start (k1), twice at provisional midpoints (k2, k3), and once at a provisional endpoint (k4), then commits the weighted sum k1/6 + k2/3 + k3/3 + k4/6. Between stages the stocks are temporarily set to the provisional states, so flows are evaluated as functions at those states — which is why flows must be pure. Randomness inside a flow under RK4 is not a stochastic model; it is four inconsistent samples averaged into one step. Put randomness in coordinate() under Euler instead (section 7).

The clock advances time to the midpoint and endpoint during the stages, so a time-dependent flow (one that reads this.currentTime) is also sampled correctly.

6. Worked example: Lotka-Volterra with auxiliaries

The shipped lotka app implements the predator-prey system

dPrey/dt = r * Prey          - g * Prey * Pred
dPred/dt = c * g * Prey * Pred - m * Pred

with visitMode: "SDRK4", dt: 0.5, maxTime: 300 in its project. The model core:

init() {
  super.init();
  this.pred = Stock.create(this.IPD, this);
  this.prey = Stock.create(this.IPY, this);
  this.prey.flow = () => this.births() - this.predation();
  this.pred.flow = () => this.eating() - this.death();
}

setup() {
  super.setup();
  this.pred.initialize(this.IPD);
  this.prey.initialize(this.IPY);
  this.showAll();
}

update() {
  super.update();
  this.showAll();
}

births()    {return this.ARE * this.prey.get();}
predation() {return this.GEE * this.prey.get() * this.pred.get();}
eating()    {return this.CEE * this.predation();}
death()     {return this.EMM * this.pred.get();}

showAll() {
  this.show("POP", this.prey.get(), this.pred.get());
  this.show("PHS", this.prey.get(), this.pred.get());
  this.show("PDC", this.pred.get());
  this.show("PYC", this.prey.get());
}

Read it against the diagram:

• Each named method is one auxiliary. predation() is computed once conceptually and consumed by both flows — prey lose it, predators gain a converted fraction eating() = c * predation(). The feedback loops of the diagram are literally the call graph of these methods.
• setup() calls stock.initialize(param) so that after a Reset the stocks restart from the current slider values, not the values captured at launch.
• One update() feeds several displays: a time-series chart (POP), a phase plot of prey versus predator (PHS), and two readouts. Display codes come from the Project.

Expected behavior: closed orbits — populations cycle, and the PHS phase plot traces a closed loop. Euler at the same dt spirals outward (the orbit artificially gains amplitude); that visible artifact is why this app runs RK4.

7. Discrete and stochastic models

Not every compartment model is an ODE. The Studio model library's StocSir implements a discrete-time stochastic SIR: per step, each susceptible independently becomes infected with probability p_SI, each infected recovers with probability p_IR. The pattern changes three things:

init() {
  super.init();
  this.S = Store.create(() => this.N - this.I0, this);
  this.I = Store.create(() => this.I0, this);
  this.R = Store.create(0, this);
}

coordinate() {
  super.coordinate();
  const lambda = this.Beta * this.I.get() / this.N;
  this.p_SI = 1 - Math.exp(-lambda);
  this.p_IR = 1 - Math.exp(-this.Sigma);
  const new_infections = APIRoot.binomial(this.S.get(), this.p_SI),
        new_recoveries = APIRoot.binomial(this.I.get(), this.p_IR);
  this.S.next = this.S.get() - new_infections;
  this.I.next = this.I.get() + new_infections - new_recoveries;
  this.R.next = this.R.get() + new_recoveries;
}

• Store instead of Stock. A plain Store is a two-phase value with no integrator: whatever you assign to next replaces the current value at commitment. There is no dt scaling and no flow function.
• coordinate() instead of flows. coordinate() is called exactly once per tick, before commitment, and only under Euler (visitMode: "SD") — the right home for random draws and if/then logic. Assigning store.next marks the store for the same synchronized commitment stocks get, so the three assignments above still read a consistent pre-step state.
• Random variates come from the inherited APIRoot statics — binomial(n, p) here; random, normal, poisson, bernoulli, choose, multinomial, and others are listed in the SDMain reference.

The transition probabilities 1 - exp(-rate) are the standard exact discretization of a constant-hazard process over one unit step, which keeps the stochastic model's mean close to the deterministic SIR of the tutorial.

Mixed models are possible: integrate continuous stocks with flows while coordinate() (Euler) or step() applies discrete events to other stores. Keep the two kinds of state in the two kinds of container.

8. Choosing Stock, Store, or sequence stock

A Stock with no flow assigned keeps its value (its pending change is zero). A Store with no assignment keeps its value as well. Details of all three behaviors, including initialization and reset, are in the Stock and Store reference.