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

-->

Using Partial-Differential Equations in MML

Introduction

This document describes how to use partial differential equations (PDEs) within MML models. JSim currently supports only 1st degree 1 dimensional PDEs.

Prerequisites:

Contents:

First Example

JSim can solve a large class of 1st degree 1-dimensional partial differential equations (PDEs). For example:

// 1st PDE example
math pde1 {
  realDomain t, x;
  t.min=0; t.max=3; t.delta=0.5;
  x.min=0; x.max=1; x.delta=0.1;
  real u(x,t);             // PDE state variable
  when (t=t.min) u = x^2;  // initial condition
  when (x=x.min) u = t^2;  // Left hand boundary conditon (LHBC)
  when (x=x.max) u:x = 0;  // Right hand boundary condition (RHBC)
  u:t = u:x:x - u;         // state equation
}
(Java plugin required)

A properly constrained PDE variable u requires the following:

  1. u is a function of at least two domains, which we will call t and x here;
  2. u is constrained for (t=t.min), this constraint is called the initial condition;
  3. u and/or u:x are constrained at both x=x.min (LHBC) and x=x.max (RHBC). The constraints must be of the form: 
          f(u, u:x) = 0

where the function f may depend on any other fully constrained variables;

  1. u possesses an entire state equation constraint that may be factored to the form: 
          A*u:t + B*u:x:x + C = 0

where A, B and C are functions of u, u:x and/or any other fully constrained variables. Note that the state equation must include both u:t and at least one of u:x or u:x:x, or else it will be treated as an ODE. An additional requirement is that B be symbolically differentiable in x. This implies B depends only on variables with explicit algebraic constraints.

Multiple Variable Problems

As with ODEs, codependent PDEs may be solved simultaneously:

// 2 codependent PDEs
math p1 {
  realDomain t, x;
  t.min=0; t.max=3; t.delta=0.25;
  x.min=0; x.max=1; x.delta=0.1;
  real u(t,x), v(t,x);
  when (t=t.min) { u=0; v=1; }
  when (x=x.min) { u=t; v = 2*t; }
  when (x=x.max) { u=t^2; v:x = 0; }
  u:t = u:x:x + u - v;
  v:t = v:x:x + v - u;
}
(Java plugin required)

JSim currently supports 3 PDE solver numeric methods. The Toms690 and Toms731 solvers will handle the full range of PDE problems described above. The BTex PDE solver handles a more restricted set of problems. The Planner will make the BTex solver available to the JSim run-time environment only if is appropriate. In cases where it is appropriate, it usually runs much faster than the Toms methods.

IC/BC Consistency in PDEs

JSim requires consistency between PDE state variable ICs and BCs. The following example is potentially inconsistent at u(t.min, x.min) depending on the value assigned to uin(t.min) at run-time:

// IC/BC consistency problem
math icbc1 {
  realDomain t, x;
  t.min=0; t.max=3; t.delta=0.25;
  x.min=0; x.max=1; x.delta=0.1;
  real u(t,x);
  extern real uin(t); 
  when (t=t.min) u = 0; 
  when (x=x.min) u = uin; 
  when (x=x.max) u:x = 0;
  u:t = u:x:x - u:x;
}
(Java plugin required)

Should uin(t.min) be assigned to a non-zero value in model icbc1 above, the model run will be rejected. The inconsistency will cause the PDE numeric methods to generate non-sensical answers.

The situation can be remedied by either modifying the IC to be consistent with the original BC, as seen below:

// IC/BC consistency via modified IC
math icbc2 {
  realDomain t, x;
  t.min=0; t.max=3; t.delta=0.25;
  x.min=0; x.max=1; x.delta=0.1;
  real u(t,x);
  extern real uin(t); 
  when (t=t.min) u = if (x=x.min) uin else 0; 
  when (x=x.min) u = uin; 
  when (x=x.max) u:x = 0;
  u:t = u:x:x - u:x;
}
(Java plugin required)

or by modifying the LHBC to be consistent with the original IC, as seen below:

// IC/BC consistency via modified BC
math icbc3 {
  realDomain t, x;
  t.min=0; t.max=3; t.delta=0.25;
  x.min=0; x.max=1; x.delta=0.1;
  real u(t,x);
  extern real uin(t); 
  when (t=t.min) u = 0; 
  when (x=x.min) u = if (t=t.min) 0 else uin; 
  when (x=x.max) u:x = 0;
  u:t = u:x:x - u:x;
}
(Java plugin required)

General Form

The current JSim Planner solves 1-dimensional PDE's with state equations of the form: 

       T*u:t + XX*u:x:x + Z = 0

and left and right hand boundary conditions (BCs) of the form: 

       when (x=x.min) f1(u) + A1*u:x = 0;
       when (x=x.max) f2(u) + A2*u:x = 0;

where

PDEs with co-dependent state variables are blocked together and solved simultaneously. The solvers operate on blocks of PDEs, not on each one individually. T, XX, Z, A1 and A2 may be functions of any u or u:x within the block, not only the u for the particular PDE.

JSim currently has some limited support for implicitly defined PDEs, but the exact extent has not yet been well defined or documented. Expanding implicit PDE definition is a continuing JSim development effort.

Available PDE Solvers

Four PDE solvers are currently available in JSim:

  1. LSFEA - Lagrangian Sliding Fluid Element Algorithm;
  2. MacCormack - 2nd order finite difference method with optional flux-corrected transport for hyperbolic and parabolic equations (JSim version 1.6.59 and above);
  3. Toms731 - adaptive moving grid method for univariate partial differential equations;
  4. Toms690 - chebyshev polynomial method for elliptic-parabolic equations.

The above order is in roughly decreasing order of speed and increasing order of generality, although this is not always strictly true. The JSim Planner determines which solvers are appropriate for a set of PDEs. The fastest amoung them becomes the default solver. The JSim user environment allows the user to select amoung available solvers.

The LSFEA Solver

The LSFEA solver accepts state equations of the form: 

      u:t = A*u:x:x + B*u:x + C

with left and right hand BCs of the following forms: 

      when (x=x.min) u = f(x,t);
      when (x=x.max) u:x = 0;

or: 

      when (x=x.min) u:x = 0;
      when (x=x.max) u = f(x,t);

or: 

      when (x=x.min) u:x = 0;
      when (x=x.max) u:x = 0;

where

Additionally:

The MacCormack Solver

The MacCormack solver accepts state equations of the form: 

       u:t = A*u:x:x + B*u:x + C

where u is the state variable and A, B and C are functions of t, x, and u. Boundary conditions must be of the form: 

      when (x=x.min) f1*u + f2*u:x = f3
      when (x=x.max) g1*u + g2*u:x = g3

where f1,2,3 and g1,2,3 are functions of t, x and u.

For further information on this method see:

The Toms731 Solver

The Toms731 solver accepts all state and boundary equations of the general form above, subject to the requirement that XX:x can be solved analytically (i.e. XX does not depend on any extern variables).

In local testing at NSR, we have found some test problems that satisfy the general form for which Toms731 does not give accurate numerical solutions. Whether this is JSim's or Toms731's problem has yet to be determined. Work continues on this topic.

Further information on Toms731 may be found at The Netlib: Toms731 .

The Toms690 Solver

The Toms690 solver was removed from JSim as of release 1.6.70. It is retained here for users of older JSim versions.

The Toms690 solver accepts all state and boundary equations of the general form above, subject to the same restrictions as Toms731.

Further information on Toms690 may be found at The Netlib: Toms690 .

[This page was last modified 03Mar08, 3:11 pm.]

Model development and archiving support at physiome.org provided by the following grants: 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.