# 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:

- Introductory Survey of MML (required)
- Introduction to the JSim GUI (recommended)

Contents:

- VF Basics
- Quick Summary of VF Basics
- Creating Delay Lines using VFs
- Data Driven Functional Relationships and VFs
- VFs and PDE Boundaries
- Comments or Questions?

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

- Compensate for VF non-linearity by choosing an appropriately fine grid spacing.
- Protect VF calls from domain value overflow (or underflow) by using "if" clauses;
- 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.

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