Skip to content

Jacobian Computation

Jacobian matrix computation utilities for numerical integration and variational equations.

Overview

Brahe provides both analytical and numerical Jacobian computation through a unified interface. Jacobian matrices describe how a function's output changes with respect to its inputs, which is essential for orbit determination, uncertainty propagation, and trajectory optimization.

DifferenceMethod 

DifferenceMethod()

Finite difference method for numerical Jacobian approximation.

Different methods trade off accuracy vs computational cost: - Forward: O(h) error, dimension+1 function evaluations - Central: O(h²) error, 2dimension function evaluations (more accurate) - Backward*: O(h) error, dimension+1 function evaluations

Example 
1
2
3
4
5
6
7
import brahe as bh

# Use central differences (most accurate)
jacobian = bh.NumericalJacobian.central(dynamics_fn)

# Or explicitly set the method
jacobian = bh.NumericalJacobian.new(dynamics_fn).with_method(bh.DifferenceMethod.CENTRAL)

Initialize instance.

BACKWARD class-attribute 

BACKWARD: Any = DifferenceMethod.Backward

Finite difference method for numerical Jacobian approximation.

Different methods trade off accuracy vs computational cost: - Forward: O(h) error, dimension+1 function evaluations - Central: O(h²) error, 2dimension function evaluations (more accurate) - Backward*: O(h) error, dimension+1 function evaluations

Example 
1
2
3
4
5
6
7
import brahe as bh

# Use central differences (most accurate)
jacobian = bh.NumericalJacobian.central(dynamics_fn)

# Or explicitly set the method
jacobian = bh.NumericalJacobian.new(dynamics_fn).with_method(bh.DifferenceMethod.CENTRAL)

CENTRAL class-attribute 

CENTRAL: Any = DifferenceMethod.Central

Finite difference method for numerical Jacobian approximation.

Different methods trade off accuracy vs computational cost: - Forward: O(h) error, dimension+1 function evaluations - Central: O(h²) error, 2dimension function evaluations (more accurate) - Backward*: O(h) error, dimension+1 function evaluations

Example 
1
2
3
4
5
6
7
import brahe as bh

# Use central differences (most accurate)
jacobian = bh.NumericalJacobian.central(dynamics_fn)

# Or explicitly set the method
jacobian = bh.NumericalJacobian.new(dynamics_fn).with_method(bh.DifferenceMethod.CENTRAL)

FORWARD class-attribute 

FORWARD: Any = DifferenceMethod.Forward

Finite difference method for numerical Jacobian approximation.

Different methods trade off accuracy vs computational cost: - Forward: O(h) error, dimension+1 function evaluations - Central: O(h²) error, 2dimension function evaluations (more accurate) - Backward*: O(h) error, dimension+1 function evaluations

Example 
1
2
3
4
5
6
7
import brahe as bh

# Use central differences (most accurate)
jacobian = bh.NumericalJacobian.central(dynamics_fn)

# Or explicitly set the method
jacobian = bh.NumericalJacobian.new(dynamics_fn).with_method(bh.DifferenceMethod.CENTRAL)

PerturbationStrategy 

PerturbationStrategy()

Strategy for computing perturbation sizes in finite differences.

The choice of perturbation size balances truncation error (wants large h) vs roundoff error (wants small h). Different strategies suit different problems.

Example 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
import brahe as bh

# Adaptive perturbation (recommended for most cases)
jacobian = bh.NumericalJacobian.new(dynamics_fn).with_adaptive(1.0, 1.0)

# Fixed absolute perturbation for all components
jacobian = bh.NumericalJacobian.new(dynamics_fn).with_fixed_offset(1e-6)

# Percentage-based perturbation
jacobian = bh.NumericalJacobian.new(dynamics_fn).with_percentage(1e-6)

Initialize instance.

adaptive builtin 

adaptive(scale_factor: float = 1.0, min_value: float = 1.0) -> PerturbationStrategy

Create an adaptive perturbation strategy.

Parameters:

Name Type Description Default
scale_factor float

Multiplier on sqrt(ε), typically 1.0

 1.0
min_value float

Minimum reference value (prevents tiny perturbations near zero)

 1.0

Returns:

Name Type Description
PerturbationStrategy PerturbationStrategy

Adaptive perturbation strategy

fixed builtin 

fixed(offset: float) -> PerturbationStrategy

Create a fixed absolute perturbation strategy.

Parameters:

Name Type Description Default
offset float

Fixed perturbation for all state components

 required

Returns:

Name Type Description
PerturbationStrategy PerturbationStrategy

Fixed perturbation strategy

percentage builtin 

percentage(percentage: float) -> PerturbationStrategy

Create a percentage-based perturbation strategy.

Parameters:

Name Type Description Default
percentage float

Percentage of state value to use as perturbation

 required

Returns:

Name Type Description
PerturbationStrategy PerturbationStrategy

Percentage perturbation strategy

NumericalJacobian 

NumericalJacobian(dynamics_fn: Any)

Numerical Jacobian provider for dynamic-sized systems using finite differences.

Computes the Jacobian numerically by perturbing the state and evaluating the dynamics. Supports forward, central, and backward finite difference methods with various perturbation strategies.

Example 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
import brahe as bh
import numpy as np

# Simple harmonic oscillator dynamics
def dynamics(t, state):
    return np.array([state[1], -state[0]])

# Default: central differences with adaptive perturbations
jacobian = bh.NumericalJacobian.new(dynamics)

# Or with custom settings:
jacobian = bh.NumericalJacobian.forward(dynamics).with_fixed_offset(1e-6)

state = np.array([1.0, 0.0])
jac_matrix = jacobian.compute(0.0, state)
print(jac_matrix)  # [[0, 1], [-1, 0]] approximately

Initialize instance.

backward builtin 

backward(dynamics_fn: callable) -> NumericalJacobian

Create with backward finite differences.

Parameters:

Name Type Description Default
dynamics_fn callable

Function with signature (t: float, state: ndarray) -> ndarray

 required

Returns:

Name Type Description
NumericalJacobian NumericalJacobian

Jacobian provider using backward differences

central builtin 

central(dynamics_fn: callable) -> NumericalJacobian

Create with central finite differences.

Parameters:

Name Type Description Default
dynamics_fn callable

Function with signature (t: float, state: ndarray) -> ndarray

 required

Returns:

Name Type Description
NumericalJacobian NumericalJacobian

Jacobian provider using central differences

compute method descriptor 

compute(t: float, state: ndarray) -> ndarray

Compute the Jacobian matrix at the given time and state.

Parameters:

Name Type Description Default
t float

Current time

 required
state ndarray

State vector at time t

 required

Returns:

Name Type Description
ndarray ndarray

Jacobian matrix ∂f/∂x (dimension × dimension)
Example 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
import brahe as bh
import numpy as np

def dynamics(t, state):
    return np.array([state[1], -state[0]])

jacobian = bh.NumericalJacobian.new(dynamics)
state = np.array([1.0, 0.5])
jac = jacobian.compute(0.0, state)
# Expected: [[0, 1], [-1, 0]] approximately

forward builtin 

forward(dynamics_fn: callable) -> NumericalJacobian

Create with forward finite differences.

Parameters:

Name Type Description Default
dynamics_fn callable

Function with signature (t: float, state: ndarray) -> ndarray

 required

Returns:

Name Type Description
NumericalJacobian NumericalJacobian

Jacobian provider using forward differences

with_adaptive method descriptor 

with_adaptive(scale_factor: float, min_value: float) -> NumericalJacobian

Set adaptive perturbation with custom parameters.

Parameters:

Name Type Description Default
scale_factor float

Multiplier on sqrt(ε), typically 1.0

 required
min_value float

Minimum reference value

 required

Returns:

Name Type Description
NumericalJacobian NumericalJacobian

Self for method chaining

with_fixed_offset method descriptor 

with_fixed_offset(offset: float) -> NumericalJacobian

Set fixed absolute perturbation for all components.

Parameters:

Name Type Description Default
offset float

Fixed perturbation size

 required

Returns:

Name Type Description
NumericalJacobian NumericalJacobian

Self for method chaining

with_method method descriptor 

with_method(method: DifferenceMethod) -> NumericalJacobian

Set the difference method.

Parameters:

Name Type Description Default
method DifferenceMethod

Finite difference method to use

 required

Returns:

Name Type Description
NumericalJacobian NumericalJacobian

Self for method chaining

with_percentage method descriptor 

with_percentage(percentage: float) -> NumericalJacobian

Set percentage-based perturbation.

Parameters:

Name Type Description Default
percentage float

Percentage of state value (e.g., 1e-6 for 0.0001%)

 required

Returns:

Name Type Description
NumericalJacobian NumericalJacobian

Self for method chaining

AnalyticJacobian 

AnalyticJacobian(jacobian_fn: Any)

Analytical Jacobian provider for dynamic-sized systems.

Uses a user-provided function that directly computes the analytical Jacobian. This is the most accurate and efficient method when the analytical Jacobian is known.

Example 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
import brahe as bh
import numpy as np

# Simple harmonic oscillator: dx/dt = v, dv/dt = -x
# Jacobian is [[0, 1], [-1, 0]]
def jacobian_fn(t, state):
    return np.array([[0.0, 1.0], [-1.0, 0.0]])

jacobian = bh.AnalyticJacobian.new(jacobian_fn)
state = np.array([1.0, 0.0])
jac = jacobian.compute(0.0, state)
print(jac)  # [[0, 1], [-1, 0]]

Initialize instance.

compute method descriptor 

compute(t: float, state: ndarray) -> ndarray

Compute the Jacobian matrix at the given time and state.

Parameters:

Name Type Description Default
t float

Current time

 required
state ndarray

State vector at time t

 required

Returns:

Name Type Description
ndarray ndarray

Jacobian matrix ∂f/∂x (dimension × dimension)
Example 
1
2
3
4
5
6
7
8
9
import brahe as bh
import numpy as np

def jacobian_fn(t, state):
    return np.array([[0.0, 1.0], [-1.0, 0.0]])

jacobian = bh.AnalyticJacobian.new(jacobian_fn)
state = np.array([1.0, 0.5])
jac = jacobian.compute(0.0, state)

See Also