DiffSL Language

The DSL is designed to be an easy and flexible language for specifying ODE systems and is based on the idea that a ODE system can be specified by a set of equations of the form:

$$ M(t) \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = F(\mathbf{u}, t) $$

where \( \mathbf{u} \) is the vector of state variables and \( t \) is the time. The DSL allows the user to specify the state vector \( \mathbf{u} \) and the RHS function \( F \). Optionally, the user can also define the derivative of the state vector \( \mathrm{d}\mathbf{u}/\mathrm{d}t \) and the mass matrix \( M \) as a function of \( \mathrm{d}\mathbf{u}/\mathrm{d}t \) (note that this function should be linear!).

The user is also free to define an an arbitrary number of intermediate scalars and vectors of the users that are required to calculate \( F \) and \( M \).

A Simple Example

To illustrate the language, consider the following simple example of a logistic growth model:

$$ \frac{\mathrm{d}N}{\mathrm{d}t} = r N (1 - N/K) $$

where \( N \) is the population, \( r \) is the growth rate, and \( K \) is the carrying capacity.

To specify this model in DiffSL, we can write:

in = [r, k]
u_i {
  N = 0.0
}
F_i {
  r * N * (1 - N/k)
}
out_i {
  N
}

Here, we define the input parameters for our model as a vector in with the growth rate r and the carrying capacity k. We then define the state vector u_i with the population N initialized to 0.0. Next, we define the RHS function F_i as the logistic growth equation. Finally, we define the output vector out_i with the population N.

Defining tensor variables

The DiffSL language only has a single type of variable, which is a n-dimensional tensor filled with double precision floating point numbers. These tensors can be dense, sparse or diagonal, and the compiler will try to choose the representation that is most efficient, preferring diagonal and then sparse matrices over dense.

The simplest tensor is a 0th dimensional scalar. For example, to define a scalar variable \( k \) with value 1.0, we write:

k { 1.0 }

Here k is the label for the tensor, which we can later use to refer to it in expressions. In the curly brackets we have one or more elements of the tensor, and here k only has a single element, which is a constant value 1.

Lets now define a 1-dimensional vector variable \( \mathbf{v} \) with 3 elements:

v_i {
 1.0,
 2.0,
 3.0,
}

The list of elements within a tensor are deliniated with a comma , and the trailing comma at the end of the list is optional. Whitespace is ignored so you can also write this tensor all on the one line:

v_i { 1.0, 2.0, 3.0 }

Subscripts

In the previous vector v_i, the subscript _i indicates that this is a 1D vector. Each subscript is a single character, and the number of subscripts indicates the number of dimensions of the tensor. You can use any character for the subscript,

v_x { 1.0, 2.0, 3.0 }
w_a { 2.0, 3.0 }
v_i { 1.0, 2.0, 3.0, 4.0 }

Ranges

Each element of a tensor can optionally give a range of index numbers, which is used by the compiler to determine the extent of each element. This is useful when defining higher dimensional tensors, such as matrices. For example, to define a 2x3 matrix \( A \) with all elements set to 1.0, we write:

A_ij {
 (0:2, 0:3) = 1.0,
}

Note the two subscript to indicate that this is a 2D tensor. The size of the single element is given in the brackets, we have two ranges 0:2 and 0:3 that correspond to the two dimensions of the matrix.

Here is another example of a 4x2 matrix \( B \) with rows 0 to 2 set to 1.0 and rows 3 to 4 set to 2.0:

A_ij {
 (0:2, 0:3) = 1.0,
 (3:4, 0:3) = 2.0,
}

For specifying a single index, you can simply write the index number without the colon, for example to define a 3x3 identity matrix $I$:

I_ij {
 (0, 0) = 1.0,
 (1, 1) = 1.0,
 (2, 2) = 1.0,
}

Note that the compiler will automatically infer the size of the tensor from the ranges you provide, so you don't need to specify the size of the tensor explicitly. Since the maximum index in the range is 2, the compiler will infer that the size of the tensor is 3x3.

Notice also that we have not defined all the elements of the matrix, only the non-zero elements. The compiler will assume that all other elements are zero.

Finally, you can also use the .. operator to specify a diagonal range of indices. For example, to define a 3x3 matrix \( D \) with the diagonal elements set to 1.0:

D_ij {
 (0..2, 0..2) = 1.0,
}

Sparse and diagonal matrices

We can automatically define a sparse matrix \( B \) by simply specifying the non-zero elements:

B_ij {
 (0, 0) = 1.0,
 (0, 1) = 2.0,
 (1, 1) = 3.0,
}

The compiler will infer that this is a 2x2 matrix, and will automatically represent it as a sparse matrix. We can force the compiler to use a dense representation by specifying the zeros explicitly:

B_ij {
 (0, 0) = 1.0,
 (0, 1) = 2.0,
 (1, 0) = 0.0,
 (1, 1) = 3.0,
}

As well as specifying a sparse matrix, we can also define a diagonal matrix by specifying the diagonal elements:

D_ij {
 (0, 0) = 1.0,
 (1, 1) = 2.0,
 (2, 2) = 3.0,
}

The compiler will infer that this is a 3x3 matrix, and will automatically represent it as a diagonal matrix.

Labels

Each element of a tensor can optionally be given a name or label that can be used to refer to the element in expressions.

For example, to define a vector with two named elements:

v_i {
 x = 1.0,
 y = 2.0,
}

Here we have defined a single tensor v_i with two named elements x and y. We can then refer to the individual elements in expressions, where they will act as if they were separate variables:

v_i { x = 1.0, y = 2.0 }
w_i { 2 * y, 3 * x }

Tensor Operations

We can use standard algebraic operations on variables. To refer to previously defined variables, we use the variable name, making sure to use the correct subscript if it is a vector or tensor.

For example, to define a scalar variable \( a \) as the sum of two other scalar variables \( b \) and \( c \), we write:

b { 1.0 }
c { 2.0 }
a { b + c }

The scalar a will therefore be equal to 3.0.

To define a vector variable \( \mathbf{v} \) as the sum of two other vector variables \( \mathbf{u} \) and \( \mathbf{w} \), we write:

u_i { 1.0, 2.0 }
w_i { 3.0, 4.0 }
v_i { u_i + w_i }

Notice that the index of the vectors within the expression must match the index of the output vector v_i. So if we defined v_a instead of v_i, the expression would be:

v_a { u_a + w_a }

Translations

For higher-dimensional tensors, the order of the indices in the expression relative to the output tensor is important.

For example, we can use the indices to define a translation of a matrix. Here we define a new matrix \( C \) that is the sum of \( A \) and \( B^T \)$, where \( B^T \) is the transpose of \( B \)

C_ij { A_ij + B_ji }

Notice that the indices of \( B^T \) are reversed in the expression compared to the output tensor \( C \), indicating that we are indexing the rows of \( B \) with j and the columns with i when we calculate the (i, j) element of \( C \).

Broadcasting

Broadcasting is supported in the language, so you can perform element-wise operations on tensors of different shapes. For example, the following will define a new vector \( \mathbf{d} \) that is the sum of \( \mathbf{a} \) and a scalar \( k \):

a_i { 1.0, 2.0 }
k { 3.0 }
d_i { a_i + k }

Here the scalar \( k \) is broadcast to the same shape as \( \mathbf{a} \) before the addition. The output vector \( \mathbf{d} \) will be \( [4.0, 5.0] \).

DiffSL uses the same broadcasting rules as NumPy, and you can read more about this in the NumPy documentation.

Contractions

The DiffSL indexing notation allows for tensor contractions, which are operations that sum over one or more indices. The rule used is that any indices that do not appear in the output tensor will be summed over.

For example, the following defines a new vector \( \mathbf{v} \) that is the sum of the rows of a matrix \( A \):

v_i { A_ij }

Another example, the following will define a new scalar \( s \) that is the sum of the element-wise product of two vectors \( \mathbf{u} \) and \( \mathbf{v} \) (i.e. an inner product):

u_i { 1.0, 2.0 }
v_i { 3.0, 4.0 }
s { u_i * v_i }

Here the i index is summed over, so the scalar s is the sum of the element-wise product of the two vectors u and v.

We can also define a matrix-vector multiplication, the following will define a new vector \( \mathbf{v} \) that is the result of a matrix-vector multiplication of a matrix \( A \) and a vector \( \mathbf{u} \):

v_i { A_ij * u_j }

This operation is actually a combination of a broadcast A_ij * u_j, followed by a contraction over the j index, the A_ij * u_j expression broadcasts the vector u to the same shape as A, forming a new 2D tensor, and the output vector v_i implicitly sums over the missing j index to form the final output vector. To illustrate this further, lets manually break this matrix-vector multiplication into two steps using an intermediary tensor M_ij:

M_ij { A_ij * u_j }
v_i { M_ij }

The first step calculates the element-wise product of A and u using broadcasting into the 2D tensor M, and the second step uses a contraction to sum over the j index to form the output vector v.

Defining a system of ODEs

The primary goal of the DiffSL language is to define a system of ODEs in the following form:

$$ \begin{align*} M(t) \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} &= F(\mathbf{u}, t) \\ \mathbf{u}(0) &= \mathbf{u}_0 \end{align*} $$

where \( \mathbf{u} \) is the vector of state variables, \( \mathbf{u}_0 \) is the initial condition, \( F \) is the RHS function, and \( M \) is the mass matrix. The DSL allows the user to specify the state vector \( \mathbf{u} \) and the RHS function \( F \).

Optionally, the user can also define the derivative of the state vector \( \mathrm{d}\mathbf{u}/\mathrm{d}t \) and the mass matrix \( M \) as a function of \( \mathrm{d}\mathbf{u}/\mathrm{d}t \) (note that this function should be linear!). The user is also free to define an arbitrary number of intermediate tensors that are required to calculate \( F \) and \( M \).

Defining state variables

To define the state variables \(\mathbf{u} \), we create a special vector variable u_i. Note the particular name u is used to indicate that this is the state vector.

The values that we use for u_i are the initial values of the state variables at \( t=0 \), so an initial condition \( \mathbf{u}(t=0) = [x(0), y(0), z(0)] = [1, 0, 0] \) is defined as:

u_i {
  x = 1,
  y = 0,
  z = 0,
}

Since we will often use the individual elements of the state vector in the RHS function, it is useful to define them as separate variables as well.

We can optionally define the time derivatives of the state variables, \( \mathbf{\dot{u}} \) as well:

dudt_i {
  dxdt = 1,
  dydt = 0,
  dzdt = 0,
}

Here the initial values of the time derivatives are given. In many cases any values can be given here as the time derivatives of the state variables are calculated from the RHS. However, if there are any algebraic variables in the equations then these values can be used used as a starting point to calculate a set of consistent initial values for the state variables.

Note that there is no need to define dudt if you do not define a mass matrix \( M \).

Defining the ODE system equations

We now define the right-hand-side function \( F \) that we want to solve, using the variables that we have defined earlier. We do this by defining a vector variable F_i that corresponds to the RHS of the equations.

For example, to define a simple system of ODEs:

$$ \begin{align*} \frac{dx}{dt} &= y \\ \frac{dy}{dt} &= -x \\ x(0) &= 1 \\ y(0) &= 0 \\ \end{align*} $$

We write:

u_i { x = 1, y = 0 }
F_i { y, -x }

Defining the mass matrix

We can also define a mass matrix \( M \) by defining a vector variable M_i which is the product of the mass matrix with the time derivative of the state vector \( M \mathbf{\dot{u}} \). This is optional, and if not defined, the mass matrix is assumed to be the identity matrix.

Notice that we are defining a vector variable M_i, which is the LHS of the ODE equations \( M \mathbf{\dot{u}} \), and not the mass matrix itself.

For example, lets define a simple DAE system using a singular mass matrix with a zero on the diagonal:

$$ \begin{align*} \frac{\mathrm{d}x}{\mathrm{d}t} &= x \\ 0 &= y-x \\ x(0) &= 1 \\ y(0) &= 0 \\ \end{align*} $$

We write:

u_i {
 x = 1,
 y = 0,
}
dudt_i {
 dxdt = 0,
 dydt = 1,
}
M_i {
 dxdt,
 0,
}
F_i {
 x,
 y-x,
}

Inputs & Outputs

Often it is useful to parameterize the system of equations using a set of input parameters. It is also useful to be able to extract certain variables from the system for further analysis. In this section we will show how to specify inputs and outputs in the DiffSL language.

Specifying inputs

We can override the values of any scalar variables by specifying them as input variables. To do this, we add a line at the top of the code to specify that these are input variables:

in = [k]
k { 1.0 }
u { 0.1 }
F { k * u }

Here we have specified a single input parameter k that is used in the RHS function F. The value of k is set to 1.0 in the code, but this value is only a default, and can be overridden by passing in a value at solve time.

We can use input parameters anywhere in the code, including in the definition of other input parameters.

in = [k]
k { 1.0 }
g { 2 * k }
F { g * u }

or in the intial conditions of the state variables:

in = [k]
k { 1.0 }
u_i {
  x = k,
}
F { u }

Specifying outputs

We can also specify the outputs of the system. These might be the state variables themselves, or they might be other variables that are calculated from the state variables.

Here is an example where we simply output the elements of the state vector:

u_i {
  x = 1.0,
  y = 2.0,
  z = 3.0,
}
out_i { x, y, z }

or we can derive additional outputs from the state variables:

u_i {
  x = 1.0,
  y = 2.0,
  z = 3.0,
}
out { x + y + z }

Pre-defined functions

The DiffSL supports the following mathematical functions that can be used in an expression:

• pow(x, y) - x raised to the power of y
• sin(x) - sine of x
• cos(x) - cosine of x
• tan(x) - tangent of x
• exp(x) - exponential of x
• log(x) - natural logarithm of x
• sqrt(x) - square root of x
• abs(x) - absolute value of x
• sigmoid(x) - sigmoid function of x
• heaviside(x) - Heaviside step function of x

You can use these functions as part of an expression in the DSL. For example, to define a variable a that is the sine of another variable b, you can write:

b { 1.0 }
a { sin(b) }

Pre-defined variables

The only predefined variable is the scalar t which is the current time, this allows the equations to be written as functions of time. For example

F_i {
  t + sin(t)
}