This page will look better in a graphical browser that supports web standards, but is accessible to any browser or internet device.

Served by Samwise.

Using Variable Functions in MML

This page is for the current JSim version 2.0. Click here for the earlier JSim 1.6 version.

This document explains variable functions (VFs), which are MML's syntax for treating a model variable v as a function v(). For example, values of a dynamic scalar variable v(t) may queried via v(expression), such as v(t-delay) or v(w(t)). This document discusses VFs in their full generality. Those whose interest is limited to creating delay lines, should instead see the document Using Delay Lines in MML.

Prerequisites:

Contents:

VF Basics

In MML, declaring a variable 1-dimensional v implies a function v(expr), where "expr" is an algebraic expression of other model variables, that can be used in model equations. Similarly, an N-dimensional variable v implies a function v(expr1, expr2, ... exprN). In the following example, w is calculated using the VF call v(t^2):

// Variable Function introductory example
math vfb1 {
  realDomain t;
  t.min=0; t.max=1; t.delta=0.1;
  real u(t) = sin(10*t);
  real v(t) = u(t^2);
}
(Java plugin required)

When the model runs, each value for v(t) will be calculated by linearly interpolating u values at the t grid points closest to t^2. Since sin(t) is not linear, the interpolation introduces some error. v(t) is a better approximation of sin(t^2) when t.delta is small than when t.delta is large. In general, VFs of non-linear expressions always introduce some interpolation error, which may be minimized by choosing an appropriate t.delta.

Using the default time grid, the required interpolation values are always well-defined because 0 < t < 1 implies t^2 is also. However, if we set t.max=2, running the model will generate an error because the required interpolation value u(4) is out of bounds. We can prevent this run-time error by using an "if" clause:

// Variable Function with grid-range protection
math vfb1 {
  realDomain t;
  t.min=0; t.max=1; t.delta=0.1;
  real u(t) = sin(10*t);
  real v(t) = if (t>1) u(1) else u(t^2);
}
(Java plugin required)

Another type of problem occurs if you set t.min=-1, all the values of v(t) for t < 0 are NaNs. This is because the calculations of u and v proceed in tandem from t.min to t.max. At t=t.min=-1, calculating v(-1) requires u(1), which hasn't been calculated yet, thus resulting in a NaN. In an ideal world, the JSim compiler would recognize this calculation order dependency and reorder the run-time calculations to accomodate it. In reality, unfortunately, this sort of order recognition is a very hard problem in complex models, and is presently without satisfactory solution. For now, modeler must be aware of calculation order dependency when writing VFs. One simple way to do this is by following the examples in this document.

Quick Summary of VF Basics

  1. Compensate for VF non-linearity by choosing an appropriately fine grid spacing.
  2. Protect VF calls from domain value overflow (or underflow) by using "if" clauses;
  3. Be aware of calculation order dependencies, following examples in this document whenever possible.

Creating Delay Lines using VFs

VFs are one of several mechanism for creating delay lines. All are discussed in the document Using Delay Lines in MML.

Data Driven Functional Relationships and VFs

VFs also provide a mechanism for incorporating external (often experimental) data into a model. Consider an initially charged capacitor(C), draining its charge(Q) over a coupled resistor(R). A simple model for this system, assuming constant R, is as follows:

// RC circuit, R constant
import nsrunit;
unit conversion on;
math rc1 {
  realDomain t sec;     // time
  t.min=0; t.max=5; t.delta=0.1;
  real R = 1 ohm;       // resistance
  real C = 1 farad;     // capacitance
  real Q0 = 1 coulomb;  // initial charge on capacitor
  real Q(t) coulomb;    // charge on capacitor at time t
  real V(t) = Q/C;      // voltage drop over capacitor
  when (t=t.min) Q=Q0;
  Q:t = -V/R;           // charge dissipated over resistor
}

Now suppose R is a non-linear function of the applied voltage that has been measured experimentally and stored in a JSim-compatible data set (see Data Files and Project Data Sets). MML for this system is below:

// RC circuit, non-linear R driven by external data
import nsrunit;
unit conversion on;
math rc1 {
  realDomain t sec;     // time
  t.min=0; t.max=5; t.delta=0.1;  // time grid definition
  real C = 1 farad;     // capacitance
  real Q0 = 1 coulomb;  // initial charge on capacitor
  realDomain vr volt;   // voltage over resistor
  vr.min=0; vr.max=Q0/C; vr.ct=101; // vr grid defintion
  extern real R(vr) ohm;// resistance
  real Q(t) coulomb;    // charge on capacitor at time t
  real V(t) = Q/C;      // voltage drop over capacitor
  when (t=t.min) Q=Q0; 
  Q:t = -V/R(V);        // charge dissipated over resistor
}
(Java plugin required)

This model introduces a new realDomain vr representing the voltage over R and that R itself is now declared "extern" with domain vr. This construct allows R's values to be drawn from the external data set via a JSim function generator. Furthermore, the ODE for Q uses the VF R(V) to calculate the resistance appropriate to the current voltage. Also note that the current voltage V(t) is a "real", not a "realDomain", and so R cannot declare V as its domain.

VFs and PDE Boundaries

Models containing PDEs in series can exhibit a VF-related problem with NaNs if boundary conditions are not properly handled. Consider flow through a stirred tank C2(t) and an spatial distributed pipe C1(t,x). The pipe's input (C1in) is drawn from the tank, and the pipe's output (C1out) empties back into the tank. A first attempt at coding might look as follows:

// tank/pipe system with pipe output calculated using VF
import nsrunit;
unit conversion on;
math vfbc1 {
  realDomain t sec;             // time grid
  t.min=0; t.max=10; t.delta=.1;  
  realDomain x cm;              // pipe spatial grid
  real L=1 cm; x.min=-L; x.max=L; x.ct=51; 
  real C1(t,x) mM, C2(t) mM;    // pipe (C1) & tank(C2) conc.
  real C1in(t) mM, C1out(t) mM; // pipe input & output conc.
  real F = 1 ml/sec;            // flow rate
  real V1 = 1 ml, V2 = 1 ml;    // pipe & tank volumes  
  real D = 1 cm^2/sec;          // diffusion constant
  when (t=t.min) { C1=C2*(x/x.max)^2; C2=1; }// initial pipe & tank conc.
  when (x=x.min) C1=C1in;       // pipe left-hand BC
  when (x=x.max) C1:x = 0;      // pipe right-hand BC
  C1:t = D*C1:x:x - 2*F*L/V1*C1:x;// pipe PDE
  C2:t = F/V2*(C1out-C2);       // tank ODE
  C1in = C2;                    // pipe input calculate from C2
  C1out = C1(t, x.max);         // pipe output calculated via VF
}

Here we calculate C1out at the pipe's right-hand boundary using a VF C1out=C1(t,x.max). Unfortunately, this code will fail at run-time with a "NaNs detected" message because the VF doesn't provide enough information to the compiler to sort out the circular dependencies of C1 and C2. If, however, we use a "when" clause to calculate C1out, the compiler can recognize the C1/C2 dependencies, and order the calculations properly:

// tank/pipe system with pipe output calculated using "when" clause
import nsrunit;
unit conversion on;
math vfbc2 {
  realDomain t sec;             // time grid
  t.min=0; t.max=10; t.delta=.1;  
  realDomain x cm;              // pipe spatial grid
  real L=1 cm; x.min=-L; x.max=L; x.ct=51; 
  real C1(t,x) mM, C2(t) mM;    // pipe (C1) & tank(C2) conc.
  real C1in(t) mM, C1out(t) mM; // pipe input & output conc.
  real F = 1 ml/sec;            // flow rate
  real V1 = 1 ml, V2 = 1 ml;    // pipe & tank volumes  
  real D = 1 cm^2/sec;          // diffusion constant
  when (t=t.min) { C1=C2*(x/x.max)^2; C2=1; }// initial pipe & tank conc.
  when (x=x.min) C1=C1in;       // pipe left-hand BC
  when (x=x.max) C1:x = 0;      // pipe right-hand BC
  C1:t = D*C1:x:x - 2*F*L/V1*C1:x;// pipe PDE
  C2:t = F/V2*(C1out-C2);       // tank ODE
  C1in = C2;                    // pipe input calculate from C2
  when (x=x.max) C1out = C1;    // pipe output calculated via "when" clause
}

This model will run without the NaN error. Here, the proper solution is to not use a VF, but to use a "when" clause instead.

Comments or Questions?

[This page was last modified 06Jul12, 3:13 pm.]

Model development and archiving support at physiome.org provided by the following grants: NIH U01HL122199 Analyzing the Cardiac Power Grid, 09/15/2015 - 05/31/2020, NIH/NIBIB BE08407 Software Integration, JSim and SBW 6/1/09-5/31/13; NIH/NHLBI T15 HL88516-01 Modeling for Heart, Lung and Blood: From Cell to Organ, 4/1/07-3/31/11; NSF BES-0506477 Adaptive Multi-Scale Model Simulation, 8/15/05-7/31/08; NIH/NHLBI R01 HL073598 Core 3: 3D Imaging and Computer Modeling of the Respiratory Tract, 9/1/04-8/31/09; as well as prior support from NIH/NCRR P41 RR01243 Simulation Resource in Circulatory Mass Transport and Exchange, 12/1/1980-11/30/01 and NIH/NIBIB R01 EB001973 JSim: A Simulation Analysis Platform, 3/1/02-2/28/07.