Jacobians and Symbolic Differentiation

Implicit algorithms and Rosenbrock-W methods need derivatives of the right-hand-side function \(f(t, x)\). This page describes how CuBIE obtains them.

What is the Jacobian?

The Jacobian matrix \(J\) of an \(n\)-state ODE system has entries:

\[J_{ij} = \frac{\partial f_i}{\partial x_j}\]

Element \(J_{ij}\) tells us how the rate of change of state \(i\) responds to a perturbation in state \(j\). Implicit solvers need the product \(J\, v\) (a Jacobian–vector product, or JVP) inside their Krylov iterations (see Nonlinear Solvers).

Some solvers compute \(J\) numerically by finite differences. This is fragile for stiff systems and expensive when \(n\) is large. CuBIE instead uses symbolic differentiation.

Symbolic Differentiation Pipeline

When you create a SymbolicODE, CuBIE:

  1. Parses the equations into SymPy expressions.

  2. Differentiates each \(f_i\) with respect to every state \(x_j\) symbolically.

  3. Applies algorithmic chain-rule steps to group common subexpressions, producing an efficient JVP function \((x, v) \mapsto J\,v\) rather than an explicit matrix.

  4. Generates CUDA-compatible code strings and writes them to a cached file.

The result is an exact, matrix-free JVP that runs as compiled device code.

Auxiliary Caching

Many subexpressions are shared between the JVP evaluation and the right-hand-side evaluation. CuBIE identifies these shared terms and caches them in a prepare step ("prepare_jac" solver helper). The subsequent JVP call ("calculate_cached_jvp") reuses the cached intermediates, avoiding redundant work.

The solver helpers are retrieved through get_solver_helper():

"linear_operator"

Basic \(J\,v\) product (no caching).

"linear_operator_cached"

\(J\,v\) product that reuses prepared intermediates.

"prepare_jac"

Evaluates and stores shared subexpressions.

"calculate_cached_jvp"

Computes \(J\,v\) using the prepared cache.

Time Derivative for Rosenbrock-W

Rosenbrock-W methods additionally require the explicit time derivative of \(f\):

\[\frac{\partial f}{\partial t} + \sum_i \frac{\partial f}{\partial d_i}\, \dot{d}_i\]

where the \(d_i\) are time-dependent driver terms. CuBIE generates this as the "time_derivative_rhs" solver helper.