v0.8.3 tutorial · numerical systems

Differentiate roots, integrals, ODEs, and DAEs

Use typed Java callbacks directly or pass a Stan function name to the supported modeling-language bindings. Parameter sensitivities are attached to the script's reusable reverse tape.

Java API + script callbacks

Roots and the implicit-function theorem

Solve F(x, θ) = 0 and propagate dx/dθ

SensitivityResult root = AlgebraicSolver.solveWithSensitivities(
    (x, theta, data, residual) ->
        residual[0] = x[0] * x[0] - theta[0],
    new double[] {1}, new double[] {2}, null,
    AlgebraicSolver.Options.defaults());
double sqrtTwo = root.values()[0][0];
double derivative = root.sensitivities()[0][0][0];

The solver uses damped Newton steps and computes dx/dθ = -Jx⁻¹Jθ. Script models call algebra_solver or algebra_solver_newton with the checked legacy signature (function, initial, theta, x_r, x_i); the result is an atomic reverse node.

Higher-order scalar integration

Pass an integrand by name

functions {
  real kernel(real x, real xc, array[] real theta,
              array[] real x_r, array[] int x_i) {
    return x_r[1] * exp(-theta[1] * x);
  }
}
// inside a procedural block
real area = integrate_1d(kernel, 0, 1, theta, x_r, x_i, 1e-8);

JDistlib uses its QUADPACK-compatible integrator, finite-differences parameter sensitivities, and Leibniz endpoint derivatives. xc is the distance to the nearest finite endpoint and is NaN for an infinite interval. The optional seventh argument is relative tolerance.

Non-stiff and stiff dynamics

Choose RK45 or adaptive BDF1

functions {
  vector decay(real t, vector y, vector theta,
               array[] real x_r, array[] int x_i) {
    return [-theta[1] * y[1]]';
  }
}
array[2] vector[1] rk = ode_rk45(decay,[1]',0,times,theta,x_r,x_i);
array[2] vector[1] stiff = ode_bdf(decay,[1]',0,times,theta,x_r,x_i);

OdeSolver provides adaptive Dormand–Prince 5(4) and forward sensitivity equations. StiffOdeSolver provides A-stable backward-Euler/BDF1 steps with adaptive step doubling; its script path computes parameter and initial-state sensitivities before attaching the trajectory to the reverse tape. The Java APIs expose tolerances and step limits through their Options classes.

Implicit index-1 dynamics

Write a residual F(t,y,y′,θ) = 0

functions {
  vector residual(real t, vector y, vector yd, vector theta,
                  array[] real x_r, array[] int x_i) {
    return [yd[1] + theta[1] * y[1]]';
  }
}
array[2] vector[1] path = dae(
    residual,[1]',[-theta[1]]',0,times,theta,x_r,x_i);

The supplied initial derivative is checked for consistency. DaeSolver advances with implicit Euler and solves each nonlinear step with AlgebraicSolver. Parameter and initial-state sensitivities are propagated to the script result.

Holonomic index-3 mechanics

Project positions and velocities onto the constraint manifold

HigherIndexDaeSolver.Result path = HigherIndexDaeSolver.integrate(
    pendulumSystem, new double[] {1,0}, new double[] {0,1},
    0, times, .002, 1e-10, 20, parameters, data);

HolonomicDaeSystem supplies unconstrained acceleration and position constraints. Projected velocity-Verlet applies Newton position projection and tangent velocity projection. Inspect maximumConstraintError(). This Java API targets mechanical holonomic systems; it is not a general DAE index-reduction engine.

Validation and boundary

Use the algorithm that matches the equations

Time origins/grids, real/integer data arrays, solver controls, and algebraic initial guesses are checked as data-only; states and parameter arrays may depend on model parameters. Independent CSV corpora cover exponential sensitivities, a stiff tracking equation, and the index-3 circular pendulum. Run ./gradlew test --tests jdistlib.StanSolverCompletionTest and ./gradlew validateModelScripts. Current boundaries are modern variadic callback signatures, event handling, adjoint sensitivities, higher-order BDF/Rosenbrock methods, and automatic index reduction for arbitrary DAEs. See fixtures 38–41 under examples/stan and the Java companion StanSolverExamples.java.