Using Integrals and Sums in MML

Introduction

This document describes using JSim MML's built-in integral() and sum() functions.

Prerequisites:

Contents:

Integrals

Integration of an integrand "expr" over realDomain t from t=lowbound to t=hibound is written as follows:

      integral(t=lowbound to hibound, expr)

Integration over the entire range of t (from t.min to t.max) may be written compactly as:

      integral(expr@t) // same as integral(t=t.min to t.max, expr)

Simple examples of integral() are demonstrated below:

// simple examples of integral() operator
unit conversion on;
import nsrunit;
math integral1 {
  realDomain t sec; t.min=0; t.max=2; t.delta=1; // time (seconds)
  real u(t) = t^2 * (1 m/sec^2);           // (meters)
  real v(t) = integral(t=t.min to t, u);   // (meter*sec)
  real w1 = integral(t=t.min to t.max, u); // (meter*sec)
  real w2 = integral(u@t);                 // (meter*sec) 
 // Another approach (solution is ODE dependent):
  real z(t);  
  when(t=t.min) z=0;  // Initial condition
  z:t = u;            // integral of u(t) from t.min to t.max
}
 

Here v(t) represents the integral of u up to the current time. w1 and w2 are calculated identically and represent the integral u over all time. All integrated variables will be assigned units m*sec, since an integral's units are the product of the units of the integrated quantity and the units of the variable of integration. Like other MML functions, integral() will generate unit conversion errors if it is used in an inappropriate context. The integral() is a first order approximation and approaches the true integral of u as t.delta -> 0 . This can be shown by letting z(t) represent the integral u over all time, expressed as an ordinary differential equation (ODE), and running the model at ever smaller t.delta values. Note that the result of an integral obtained by an ODE equation is solver dependent (Euler accuracy is usually much worse than Radau).

While integrating from t.min to t.max or from t.min to t are the most common uses, the integrals may be between any two values to t. Also, the expression integrated need not be a simple variable, for example:

      integral(t=t.min+1 to t/2, u^2*exp(-v))

MML integral() can also be used to "integrate out" one domain of a multi-dimensional variable:

// "integrating out" one variable domain
unit conversion on;
import nsrunit;
math integral2 {
  realDomain t sec; t.min=0; t.max=6; t.delta=1; // time (seconds)
  realDomain x m; x.min=0; x.max=1; x.delta=0.1; // space (meters)
  real u(t,x) = t^2*(1 mole/sec^2) + x^2*(1 mole/m^2); // u (mole)
  real v(t) = integral(x=x.min to x.max, u); // v (mole*meters)
}


 

 

Sums

Using the sum and integral operators to integrate a function IS NOT RECOMMENDED and is illustrated by the following example.

import nsrunit;
unit conversion on;
math sum0{
realDomain t sec; t.min=0; t.max=10; t.ct=101;
real x(t) = t^2;

/* Compute the average over x(t):
   x(t) is a discrete series. */
   real Average sec^2;
   Average = sum(x@t)/t.ct;

/* An approximate integral can be calculated in
   two different ways: The first way,
   ApproximateIntegral1 uses the sum operator
   AND IS NOT RECOMMENDED;   */
   real ApproximateIntegral1 sec^3;
   ApproximateIntegral1 = sum(x*t.delta@t);
// yields ApproximateIntegral1=338.35.

/* The second way is more accurate. It uses
   the integral operator which uses trapezoidal integration.
*/
   real  ApproximateIntegral2 sec^3;
   ApproximateIntegral2 = integral(x@t);
// yields ApproximateIntegral2=333.35

/* Both of these methods converge to the exact answer as
   t.ct (the number of points between t.min=0 and t.max=10
   gets large. The most accurate method is to directly
   integrate the function.                                */
   real ExactIntegral(t) sec^3;
   when(t=t.min) ExactIntegral = 0;
   ExactIntegral:t = x;
// yields ExactIntegral = 333.33333333333.
}
 

The sum operator is used extensively in the models found by searching with the keyword "Fourier."

The summation of a summand "expr" over realDomain t from t=lowbound to t=hibound is written as follows:

      sum(t=lowbound to hibound, expr)

Summation over the entire range of t (from t.min to t.max) may be written compactly as:

      sum(expr@t) // same as sum(t=t.min to t.max, expr)

Simple examples of sum() are demonstrated below:

// simple examples of sum() operator
unit conversion on;
import nsrunit;
math sums {
  realDomain t sec; t.min=0; t.max=5; t.delta=1; // time (seconds)
  real u(t) = t^2 * (1 m/sec^2);  // u (meters)
  real v(t) = sum(t=t.min to t, u); // v (meter)
  real w1 = sum(t=t.min to t.max, u); // w1 (meter)
  real w2 = sum(u@t); // w2 (meter)
}
 

Here v(t) represents the sum of u up to the current time. w1 and w2 are calculated identically and represent the sum u over all time. All summed variables will be assigned units m, since a sum's units are equal to those of the summand. Like other MML functions, sum() will generate unit conversion errors if it is used in an inappropriate context.

While summing from t.min to t.max or from t.min to t are the most common uses, sums may be between any two values to t. Also, the summand need not be a simple variable, for example:

      sum(t=t.min+1 to t/2, u^2*exp(-v))

MML () can also be used to "sum out" one domain of a multi-dimensional variable:

// "summing out" one variable domain
unit conversion on;
import nsrunit;
math sum2 {
  realDomain t sec; t.min=0; t.max=6; t.delta=1;
  realDomain x m; x.min=0; x.max=1; x.delta=0.1;
  real u(t,x) = t^2*(1 mole/sec^2) + x^2*(1 mole/m^2); // u in mole
  real v(t) = sum(x=x.min to x.max, u); // v in mole
}
 

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.