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. Requesting any Jacobian-carrying helper with the "cached" variant returns a result that reuses the prepared intermediates and carries the prepare_jac device function that fills the cache.

The solver helpers are retrieved through get_solver_helper() with a role and a variant:

role="linear_operator"

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

role="linear_operator", jacobian_at="step"

\(J\,v\) product that reuses prepared intermediates; the returned result also carries prepare_jac, which evaluates and stores the shared subexpressions.

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.