Implicit Equations in MML

Introduction

This document describes using implicit equations within JSim MML models. Implicit equations refer to equations in which the unknown variable(s) are not isolated on one side of the equation, but are intermixed with other known and unknown variables.

Implicit equations may be linear or non-linear and may involve one or many unknown variables. The number of equations must be the same as the number of unknowns. If it is less, the problem will be rejected as underspecified. If it is more, the problem will be rejected as overspecified. (There are cases where N+1 equations exactly specify N unknowns, but JSim is not yet smart enough to recognize them.) Unknown variables may have any number of domains, however all variables that are solved together must have the same set of domains.

Prerequisites:

Contents:

  • Linear Implicit Equations in One Variable
  • Linear Implicit Equations of Multiple Variables
  • Linear Implicit Equations in ODEs
  • Non-linear Implicit Equations
  • Bounded and Unbounded Non-Linear Implicit Equations
  • The Approximately-Equals Relation
  • Non-Linear Solver Run-Time Controls
  • Non-linear Implicit Equations in Multiple Variables
  • Limitations
  • Comments or Questions?
    Give feedback

Linear Implicit Equations in One Variable

Linear implicit equations of a single variable are solved internally by JSim via algebraic variable isolation. JSim solves the following:

// linear implicit eqn in one variable
math implicit1 {
  real a=1, b=2, y=4;
  real x;
  a*x - b = y;
}
 

 

by internally rearranging it as:

 

      x = (y+b)/a;

 

Therefore, model writers should not hesitate to use the implicit form of such equations if they express the underlying ideas more clearly.

Linear Implicit Equations of Multiple Variables

These problems are often call simultaneous equations. For example:

// linear implicit eqns of 2 variables
math implicit2 {
  realDomain t; t.min=0; t.max=10; t.delta=1;
  real a=4, b=5;
  real x(t), y(t);
  a*x + b*y = t;
  b*x - y = a;
}

 

JSim solves such problems via a matrix inversion operation. If the equations are not independent, the matrix will not be invertible and an error will be generated at model run time.

Linear Implicit Equations in ODEs

Linear implicit equations can be used to define the initial conditions and/or state equations for ordinary differential equations (ODEs). An example involving initial conditions:

// linear implicit eqns in ODE ICs
math implicit3 {
  realDomain t; 
  t.min=0; t.max=10; t.delta=1;
  real u(t), v(t);
  when (t=t.min) {
    u + 2*v = 4;
    u - v = 1;
  }
  u:t = v-u;
  v:t = u-v;
}
 

 

When defining ODE state equations, the derivates of the state variables are treated as unknowns, while the state variables themselves are treated as knowns. Implicit equations allow the simultaneous calculation of several derivates:

 

// linear implicit ODE state eqns
math implicit4 {
  realDomain t;
  t.min=0; t.max=10; t.delta=1;
  real u(t), v(t);
  when (t=t.min) {
    u=2;
    v=1;
  }
  u:t + v:t = u-v^2;
  u:t - v:t = u+v^2; 
}
 

Non-linear Implicit Equations

Consider a simple example with a single unknown variable (x):

      realDomain t;
      real x(t);
      x^3 - t = 0;

This example yields x as the cube-root to t. Easy pie.

However, many non-linear implicit equations are trickier than the above example because they may have no solution or multiple solutions. Consider the following system:

      realDomain t;
      real x(t);
      x^2 = t;

In this system, x has two possible solutions for all positive t (x's positive and negative square roots). There is a unique x solution at t=0, but no real x solution for negative t. (Complex variables are contemplated in JSim's future, but are not yet available). At run time this model would generate an error if t assumes any negative values. For positive tx would be assigned to either the positive or negative square root without guarantee as to which. In fact, some x may be positive while others are negative.

This unpredictability is often unacceptable in predictive modeling situations. In the example above, for instance, the modeler may always desire the non-negative solution. In this particular case, it would be trivial to rewrite the model as:

      x = sqrt(a);  // sqrt() always returns a non-negative value

however, non-linear equations are often too complicated for this approach to work.

To combat the unpredicability of multiple solutions, models may be refined by using relational expressions to bound solution values. Supported relational operators are:

  • "<" - less than;
  • "<=" - less than or equal to;
  • "~=" - approximately equal to (clarification follows);
  • ">=" - greater than or equal to;
  • ">" - greater than.

In the current case, an ideal respecification for a non-negative solution for x would be:

      realDomain t;
      real x(t);
      x^2 = t;
      x >= 0;

Unfortunately, JSim's current numeric solvers are not quite sophisticated enough to solve this problem as specified. To understand why, we need to look a little more closely at JSim's underlying numeric solvers.

Bounded and Unbounded Non-Linear Implicit Equations

The JSim numeric engine distinguishes two cases of non-linear implicit equations. Bounded problems are those in which the model specifies both a low bound (">" or ">=") and a high bound ("<" or "<=") for each variable to be solved. Unbounded problems are those in which at least one variable to be solved is missing at least one bound. Bounded problems are, by default, solved via the simplex algorithm which uses the bounds explicitly to limit its search. Unbounded problems are solved via the GGopt algorithm which does not use bounds at all.

JSim currently does not have any algorithms to solve "semi-bounded" problems such as the above square-root problem with a low bound of "x>=0" but no upper bound. As a consequence the model will be solved via the GGopt algorithm, which ignores the low bound during its calculations, and may possibly return a negative value. Therefore, the proper formulation for square-root model is to provide an upper bound in addition to a lower bound, turning it into a fully bounded problem which allows simplex to work properly. The following reformulation should work properly for all positive t:

// bounded non-linear implicit eqn
math implicit5 {
  realDomain t;
  t.min=0; t.max=10; t.delta=1;
  real x(t);
  x^2 = t;
  x >= 0;
  x <= if (t<1) 1 else t;
}
 

 

JSim guarantees that any solution returned will satisfy all specified relations. Therefore, if ggopt is used in problems containing bounds (either because the problem is semi-bound, or because the user specifies ggopt as the bounded solver control below) the returned value may or may not pass the relations tests. If the returned value does not satify the required relation, an run-time error will result. Summary of the three variants of the square-root model:

 

  • unbounded: x may return negative or positive;
  • semi-bounded: x will return positive, or an error will result;
  • bounded: x will return positive.

The Approximately-Equals Relation

The "~=" relation is currently used by JSim to give the simplex and GGopt numeric algorithms a starting value to search. A good starting guess can significantly improve the performance of these algorithms. If there multiple solutions, the starting value contributes to which one is selected. This often means that the selected solution will be the one closest to the starting value, but this is not always the case, so modelers should not rely upon it for stable models.

In the future, JSim may incorporate algorithms that find all solutions to a set of implicit equations and use the approximately-equals relation to pick the closest one. Such algorithms are not yet available, however.

Non-Linear Solver Run-Time Controls

Like the controls that govern ODE and PDE solvers, there are user-adjustable controls that affect non-linear implicit equation solvers. (There are currently no user-adjustable controls for linear implicit equations). This set of controls will be visible in the "Solvers" page of JSim's model run-time tab when a model requires them. The controls are detailed below:

  • unbound: This controls which algorithm is used for unbounded implicit equations. Currently, the only available algorithm is "ggopt", but that may change in the future;
  • bound: This controls which algorithm is used for bounded implicit equations. The default is "simplex", but "ggopt" is also available. If ggopt is selected, a run-time error will result if the returned value does not satisfy the specified bounds.
  • errtol: This is the average root-mean-square error required for the equations to be considered solved. If the solvers cannot achieve this level of accuracy, a run-time error will result. Decreasing this number will increase the accuracy of the solution, but increase the probability of run-time error. Increasing this number will decrease the probability of run-time error, but decrease solution accuracy.
  • maxcalls: This number, multiplied by the number of simultaneous variables to be solved, gives to maximum number of function calls to try in attempting to find a solution.
  • maxiters: This number, multiplied by the number of simultaneous variables to be solved, gives the maximum number of ggopt iterations (which ggopt is used).
  • eps: GGopt relative error.
  • istep: Simplex initial step size.

Non-linear Implicit Equations in Multiple Variables

These problems are a straight-forward analogy of the single variable case. The following example finds points of intersection between a simple parabola and the unit circle:

      real x, y;
      x^2 + y^2 = 1;
      y = x^2;

Note that there are two solutions, symmetric about the Y axis.

If we were interested in only the positive X value, we could make this a bound problem via:

// 2 bounded non-linear implicit eqns
math implicit6 {
  real x, y;
  x^2 + y^2 = 1;
  y = x^2;
  x >= 0;
  x <= 1;
  y >= 0;
  y <= 1;
}
 

 

Note that upper and lower bounds for each variable are required. Otherwise, the unbounded algorithm will be used, possible generating errors as explained above.

 

Limitations

JSim's current non-linear zero finders share non-reentrant native simplex and ggopt methods, meaning that models using them cannot run the model optimizer successfully. Work is under way to correct this situation.

Comments or Questions?

Give feedback

Model development and archiving support at https://www.imagwiki.nibib.nih.gov/physiome 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.