Differentiable CSTR: Sensitivity Analysis#
This notebook provides a deep dive into sensitivity analysis using automatic differentiation.
What You’ll Learn#
∂Output/∂Input - How outlet composition changes with inlet conditions
∂Output/∂Parameters - Sensitivity to kinetic parameters (A, Ea)
∂Output/∂Operating - Sensitivity to operating conditions (V, T)
Jacobian matrices - Full input-output sensitivity
Hessian matrices - Second-order sensitivities (curvature)
Gradient-based optimization - Using sensitivities to optimize
The CSTR Model#
We consider a simple CSTR with first-order kinetics:
\[A \rightarrow B\]
Rate: \(r = k \cdot C_A\) where \(k = A \cdot \exp(-E_a/RT)\)
import jax
import jax.numpy as jnp
from jax import Array
# Enable 64-bit precision
jax.config.update("jax_enable_x64", True)
import optimistix as optx
from difflow.streams import Stream, make_stream, get_flows
from difflow.thermo import IdealThermo, SpeciesData
from difflow.units.cstr import CSTR, CSTRParams
WARNING:2026-03-01 07:48:38,916:jax._src.xla_bridge:905: Platform 'METAL' is experimental and not all JAX functionality may be correctly supported!
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
W0000 00:00:1772369318.917035 12631935 mps_client.cc:510] WARNING: JAX Apple GPU support is experimental and not all JAX functionality is correctly supported!
I0000 00:00:1772369318.926279 12631935 service.cc:145] XLA service 0xa0bec1700 initialized for platform METAL (this does not guarantee that XLA will be used). Devices:
I0000 00:00:1772369318.926289 12631935 service.cc:153] StreamExecutor device (0): Metal, <undefined>
I0000 00:00:1772369318.927283 12631935 mps_client.cc:406] Using Simple allocator.
I0000 00:00:1772369318.927290 12631935 mps_client.cc:384] XLA backend will use up to 55662313472 bytes on device 0 for SimpleAllocator.
WARNING:2026-03-01 07:48:38,927:jax._src.xla_bridge:905: Platform 'mps' is experimental and not all JAX functionality may be correctly supported!
Metal device set to: Apple M4 Pro
Setup: Species and Thermodynamics#
species_data = {
"A": SpeciesData(
name="A", MW=100.0, Cp_coeffs=(75.0, 0.0, 0.0, 0.0),
Hvap_coeffs=(35000.0, 0.38, 500.0), antoine_coeffs=(10.0, 3000.0, -50.0),
Hf=0.0,
),
"B": SpeciesData(
name="B", MW=100.0, Cp_coeffs=(75.0, 0.0, 0.0, 0.0),
Hvap_coeffs=(30000.0, 0.38, 450.0), antoine_coeffs=(10.0, 2800.0, -40.0),
Hf=-50000.0,
),
}
thermo = IdealThermo(species_data)
species_order = ["A", "B"]
stoichiometry = jnp.array([[-1.0], [+1.0]])
def rate_function(C: dict[str, Array], T: Array, params: dict) -> Array:
"""First-order reaction: A → B with Arrhenius kinetics."""
k = params["A"] * jnp.exp(-params["Ea"] / (8.314 * T))
return jnp.array([k * C["A"]])
print("Setup complete ✓")
Setup complete ✓
1. Sensitivity to Inlet Conditions (∂Output/∂Input)#
How does the outlet composition change when we vary the inlet?
Question: If we increase inlet A flow by 1 mol/s, how much more B do we produce?
def outlet_B_vs_inlet(F_A_in: Array, F_B_in: Array, T_in: Array) -> Array:
"""Compute outlet F_B as function of inlet conditions."""
params = CSTRParams(
V=jnp.array(1.0),
rate_fn=rate_function,
stoich=stoichiometry,
rate_params={"A": jnp.array(1e6), "Ea": jnp.array(50000.0)},
species_order=species_order,
dH_rxn=jnp.array([-50000.0]),
)
cstr = CSTR(params, thermo=thermo, mode="isothermal")
inlet = make_stream({"A": F_A_in, "B": F_B_in}, T=T_in, P=101325.0)
outlet, _ = cstr(inlet, T_spec=350.0)
return outlet["F_B"]
# Base case
F_A_base, F_B_base, T_base = 10.0, 0.0, 300.0
F_B_out = outlet_B_vs_inlet(
jnp.array(F_A_base), jnp.array(F_B_base), jnp.array(T_base)
)
print("Base case:")
print(f" Inlet: F_A={F_A_base} mol/s, F_B={F_B_base} mol/s, T={T_base} K")
print(f" Outlet F_B = {float(F_B_out):.4f} mol/s")
Base case:
Inlet: F_A=10.0 mol/s, F_B=0.0 mol/s, T=300.0 K
Outlet F_B = 9.9480 mol/s
# Compute gradients with respect to each inlet variable
grad_FA = jax.grad(outlet_B_vs_inlet, argnums=0)
grad_FB = jax.grad(outlet_B_vs_inlet, argnums=1)
grad_T = jax.grad(outlet_B_vs_inlet, argnums=2)
dFB_dFA_in = grad_FA(jnp.array(F_A_base), jnp.array(F_B_base), jnp.array(T_base))
dFB_dFB_in = grad_FB(jnp.array(F_A_base), jnp.array(F_B_base), jnp.array(T_base))
dFB_dT_in = grad_T(jnp.array(F_A_base), jnp.array(F_B_base), jnp.array(T_base))
print("Sensitivities (gradients):")
print(f" ∂F_B_out/∂F_A_in = {float(dFB_dFA_in):.4f}")
print(f" ∂F_B_out/∂F_B_in = {float(dFB_dFB_in):.4f}")
print(f" ∂F_B_out/∂T_in = {float(dFB_dT_in):.6f} mol/s per K")
print("\n📊 Interpretation:")
print(f" • Increasing inlet A by 1 mol/s → outlet B increases by {float(dFB_dFA_in):.2f} mol/s")
print(f" • Inlet B passes through unchanged (∂F_B/∂F_B_in ≈ 1)")
print(f" • Inlet T has minimal effect at isothermal operation")
Sensitivities (gradients):
∂F_B_out/∂F_A_in = 0.9896
∂F_B_out/∂F_B_in = 0.9948
∂F_B_out/∂T_in = 0.000000 mol/s per K
📊 Interpretation:
• Increasing inlet A by 1 mol/s → outlet B increases by 0.99 mol/s
• Inlet B passes through unchanged (∂F_B/∂F_B_in ≈ 1)
• Inlet T has minimal effect at isothermal operation
2. Sensitivity to Kinetic Parameters#
How does conversion change with the Arrhenius parameters?
This is crucial for:
Parameter estimation: Fitting A and Ea to experimental data
Uncertainty propagation: How parameter uncertainty affects predictions
def conversion_vs_kinetics(log_A: Array, Ea: Array) -> Array:
"""Compute conversion as function of kinetic parameters.
Uses log(A) instead of A for better numerical conditioning.
"""
params = CSTRParams(
V=jnp.array(1.0),
rate_fn=rate_function,
stoich=stoichiometry,
rate_params={"A": jnp.exp(log_A), "Ea": Ea},
species_order=species_order,
dH_rxn=jnp.array([-50000.0]),
)
cstr = CSTR(params, thermo=thermo, mode="isothermal")
inlet = make_stream({"A": 10.0, "B": 0.0}, T=300.0, P=101325.0)
outlet, info = cstr(inlet, T_spec=350.0)
return info["conversion"]["A"]
# Base case
log_A_base = jnp.log(1e6) # A = 1e6 /s
Ea_base = jnp.array(50000.0) # 50 kJ/mol
X = conversion_vs_kinetics(log_A_base, Ea_base)
print(f"Base case: A = 1e6 /s, Ea = 50 kJ/mol")
print(f"Conversion = {float(X)*100:.2f}%")
Base case: A = 1e6 /s, Ea = 50 kJ/mol
Conversion = 99.48%
# Compute gradients
grad_fn = jax.grad(conversion_vs_kinetics, argnums=(0, 1))
dX_dlogA, dX_dEa = grad_fn(log_A_base, Ea_base)
print("Sensitivities:")
print(f" ∂X/∂(log A) = {float(dX_dlogA):.4f}")
print(f" ∂X/∂Ea = {float(dX_dEa)*1000:.6f} per kJ/mol")
print("\n📊 Physical interpretation:")
print(f" • Doubling A (Δlog A = 0.693) increases X by {float(dX_dlogA)*0.693*100:.2f}%")
print(f" • Increasing Ea by 1 kJ/mol changes X by {float(dX_dEa)*1000*100:.2f}%")
print("\n💡 Application:")
print(" These gradients enable gradient-based fitting of A and Ea")
print(" to experimental conversion data!")
Sensitivities:
∂X/∂(log A) = 0.0052
∂X/∂Ea = -0.001777 per kJ/mol
📊 Physical interpretation:
• Doubling A (Δlog A = 0.693) increases X by 0.36%
• Increasing Ea by 1 kJ/mol changes X by -0.18%
💡 Application:
These gradients enable gradient-based fitting of A and Ea
to experimental conversion data!
3. Sensitivity to Operating Conditions (V, T)#
Which is more effective for increasing production: larger reactor or higher temperature?
def outlet_B_vs_operating(V: Array, T_reactor: Array) -> Array:
"""Compute outlet F_B as function of V and T."""
params = CSTRParams(
V=V,
rate_fn=rate_function,
stoich=stoichiometry,
rate_params={"A": jnp.array(1e6), "Ea": jnp.array(50000.0)},
species_order=species_order,
dH_rxn=jnp.array([-50000.0]),
)
cstr = CSTR(params, thermo=thermo, mode="isothermal")
inlet = make_stream({"A": 10.0, "B": 0.0}, T=300.0, P=101325.0)
outlet, _ = cstr(inlet, T_spec=T_reactor)
return outlet["F_B"]
# Base case
V_base = jnp.array(1.0)
T_base = jnp.array(350.0)
F_B = outlet_B_vs_operating(V_base, T_base)
print(f"Base case: V = 1.0 m³, T = 350 K")
print(f"Outlet F_B = {float(F_B):.4f} mol/s")
Base case: V = 1.0 m³, T = 350 K
Outlet F_B = 9.9480 mol/s
# Compute gradients
grad_V = jax.grad(outlet_B_vs_operating, argnums=0)
grad_T = jax.grad(outlet_B_vs_operating, argnums=1)
dFB_dV = grad_V(V_base, T_base)
dFB_dT = grad_T(V_base, T_base)
print("Sensitivities:")
print(f" ∂F_B/∂V = {float(dFB_dV):.4f} mol/s per m³")
print(f" ∂F_B/∂T = {float(dFB_dT):.6f} mol/s per K")
# Elasticities (normalized sensitivities)
elasticity_V = float(dFB_dV) * float(V_base) / float(F_B)
elasticity_T = float(dFB_dT) * float(T_base) / float(F_B)
print("\n📊 Elasticities (% change in F_B per % change in parameter):")
print(f" ε_V = {elasticity_V:.4f}")
print(f" ε_T = {elasticity_T:.4f}")
print("\n (A 1% increase in V or T causes this % change in F_B)")
Sensitivities:
∂F_B/∂V = 0.0517 mol/s per m³
∂F_B/∂T = 0.002538 mol/s per K
📊 Elasticities (% change in F_B per % change in parameter):
ε_V = 0.0052
ε_T = 0.0893
(A 1% increase in V or T causes this % change in F_B)
4. Jacobian Matrix: Full Input-Output Mapping#
The Jacobian matrix contains all first-order sensitivities:
\[J_{ij} = \frac{\partial y_i}{\partial x_j}\]
For a CSTR mapping \([F_{A,in}, F_{B,in}, T_{in}] \rightarrow [F_{A,out}, F_{B,out}, T_{out}]\)
def cstr_function(inputs: Array) -> Array:
"""CSTR as a vector function: [F_A_in, F_B_in, T_in] → [F_A_out, F_B_out, T_out]."""
F_A_in, F_B_in, T_in = inputs[0], inputs[1], inputs[2]
params = CSTRParams(
V=jnp.array(1.0),
rate_fn=rate_function,
stoich=stoichiometry,
rate_params={"A": jnp.array(1e6), "Ea": jnp.array(50000.0)},
species_order=species_order,
dH_rxn=jnp.array([-50000.0]),
)
cstr = CSTR(params, thermo=thermo, mode="isothermal")
inlet = make_stream({"A": F_A_in, "B": F_B_in}, T=T_in, P=101325.0)
outlet, _ = cstr(inlet, T_spec=350.0)
return jnp.array([outlet["F_A"], outlet["F_B"], outlet["T"]])
# Compute Jacobian using forward-mode AD
inputs = jnp.array([10.0, 0.0, 300.0])
jacobian = jax.jacfwd(cstr_function)(inputs)
outputs = cstr_function(inputs)
print(f"Inputs: F_A_in={inputs[0]:.1f}, F_B_in={inputs[1]:.1f}, T_in={inputs[2]:.1f}")
print(f"Outputs: F_A_out={outputs[0]:.4f}, F_B_out={outputs[1]:.4f}, T_out={outputs[2]:.1f}")
print("\nJacobian matrix J = ∂(outputs)/∂(inputs):")
print(" dF_A_in dF_B_in dT_in")
print(f"dF_A_out {jacobian[0, 0]:9.4f} {jacobian[0, 1]:9.4f} {jacobian[0, 2]:9.6f}")
print(f"dF_B_out {jacobian[1, 0]:9.4f} {jacobian[1, 1]:9.4f} {jacobian[1, 2]:9.6f}")
print(f"dT_out {jacobian[2, 0]:9.4f} {jacobian[2, 1]:9.4f} {jacobian[2, 2]:9.6f}")
Inputs: F_A_in=10.0, F_B_in=0.0, T_in=300.0
Outputs: F_A_out=0.0520, F_B_out=9.9480, T_out=350.0
Jacobian matrix J = ∂(outputs)/∂(inputs):
dF_A_in dF_B_in dT_in
dF_A_out 0.0104 0.0052 -0.000000
dF_B_out 0.9896 0.9948 0.000000
dT_out 0.0000 0.0000 0.000000
Interpreting the Jacobian#
Row i: How output i changes with each input
Column j: How input j affects each output
Applications:
Linear uncertainty propagation: \(\text{Var}(y) = J \cdot \text{Var}(x) \cdot J^T\)
Controllability analysis: Which inputs affect which outputs?
Stability analysis: Eigenvalues of Jacobian
5. Hessian Matrix: Second-Order Sensitivities#
The Hessian tells us about the curvature of the objective:
\[H_{ij} = \frac{\partial^2 f}{\partial x_i \partial x_j}\]
This is crucial for:
Understanding diminishing returns
Newton-based optimization
Second-order uncertainty propagation
def conversion_vs_VT(params: Array) -> Array:
"""Conversion as function of [V, T]."""
V, T_reactor = params[0], params[1]
cstr_params = CSTRParams(
V=V,
rate_fn=rate_function,
stoich=stoichiometry,
rate_params={"A": jnp.array(1e6), "Ea": jnp.array(50000.0)},
species_order=species_order,
dH_rxn=jnp.array([-50000.0]),
)
cstr = CSTR(cstr_params, thermo=thermo, mode="isothermal")
inlet = make_stream({"A": 10.0, "B": 0.0}, T=300.0, P=101325.0)
outlet, info = cstr(inlet, T_spec=T_reactor)
return info["conversion"]["A"]
# Compute gradient and Hessian
params = jnp.array([1.0, 350.0])
X = conversion_vs_VT(params)
grad = jax.grad(conversion_vs_VT)(params)
hessian = jax.hessian(conversion_vs_VT)(params)
print(f"At V = {params[0]:.1f} m³, T = {params[1]:.1f} K:")
print(f" Conversion X = {float(X)*100:.2f}%")
print(f"\nGradient (first derivatives):")
print(f" ∂X/∂V = {float(grad[0]):.4f}")
print(f" ∂X/∂T = {float(grad[1]):.6f}")
print(f"\nHessian (second derivatives):")
print(f" ∂²X/∂V² = {float(hessian[0, 0]):.4f}")
print(f" ∂²X/∂T² = {float(hessian[1, 1]):.10f}")
print(f" ∂²X/∂V∂T = {float(hessian[0, 1]):.6f}")
At V = 1.0 m³, T = 350.0 K:
Conversion X = 99.48%
Gradient (first derivatives):
∂X/∂V = 0.0052
∂X/∂T = 0.000254
Hessian (second derivatives):
∂²X/∂V² = -0.0103
∂²X/∂T² = -0.0000137830
∂²X/∂V∂T = -0.000251
Interpreting the Hessian#
∂²X/∂V² < 0: Diminishing returns with volume (concave) - adding more volume gives less incremental benefit
∂²X/∂T²: Curvature of temperature dependence
∂²X/∂V∂T: Cross-effect - how does the V sensitivity change with T?
6. Gradient-Based Optimization#
Now let’s use gradients to optimize the reactor design!
Problem Formulation#
We want to find the reactor volume \(V\) and temperature \(T\) that maximize profit:
\[\text{maximize} \quad X - 0.01 \cdot \text{Cost}\]
where:
\(X\) = conversion (fraction of A reacted)
\(\text{Cost} = V + 0.001 \cdot (T - 300)^2\) (capital + energy costs)
Box Constraints#
Physical and practical limits constrain our design:
Volume: \(V \in [0.1, 5.0]\) m³ (minimum practical size, maximum space available)
Temperature: \(T \in [300, 450]\) K (ambient to material limits)
Solution Approach#
Since BFGS is an unconstrained optimizer, we use a sigmoid transformation to map unbounded variables to the constrained domain. This is a common technique for handling box constraints with gradient-based methods.
def objective(params: Array, args=None) -> Array:
"""Objective to minimize (negative of profit).
We minimize the negative because optimistix.minimise finds minima.
Minimizing -profit is equivalent to maximizing profit.
"""
V, T_reactor = params[0], params[1]
cstr_params = CSTRParams(
V=V,
rate_fn=rate_function,
stoich=stoichiometry,
rate_params={"A": jnp.array(1e6), "Ea": jnp.array(50000.0)},
species_order=species_order,
dH_rxn=jnp.array([-50000.0]),
)
cstr = CSTR(cstr_params, thermo=thermo, mode="isothermal")
inlet = make_stream({"A": 10.0, "B": 0.0}, T=300.0, P=101325.0)
outlet, info = cstr(inlet, T_spec=T_reactor)
conversion = info["conversion"]["A"]
cost = V + 0.001 * (T_reactor - 300.0) ** 2
profit = conversion - 0.01 * cost
return -profit # Minimize negative profit = maximize profit
# Box constraints
V_bounds = (0.1, 5.0) # m³
T_bounds = (300.0, 450.0) # K
def transform_params(x: Array) -> Array:
"""Transform unconstrained x to bounded [V, T] using sigmoid.
sigmoid(x) maps (-∞, +∞) → (0, 1), then we scale to bounds.
This ensures the optimizer never violates constraints.
"""
sigmoid = lambda z: 1.0 / (1.0 + jnp.exp(-z))
V = V_bounds[0] + (V_bounds[1] - V_bounds[0]) * sigmoid(x[0])
T = T_bounds[0] + (T_bounds[1] - T_bounds[0]) * sigmoid(x[1])
return jnp.array([V, T])
def bounded_objective(x: Array, args=None) -> Array:
"""Objective in transformed (unconstrained) space."""
return objective(transform_params(x))
def inverse_transform(params: Array) -> Array:
"""Map bounded params back to unconstrained space (for initial guess)."""
def logit(p):
return jnp.log(p / (1.0 - p))
v_norm = (params[0] - V_bounds[0]) / (V_bounds[1] - V_bounds[0])
t_norm = (params[1] - T_bounds[0]) / (T_bounds[1] - T_bounds[0])
return jnp.array([logit(v_norm), logit(t_norm)])
# Initial guess
initial_params = jnp.array([0.5, 320.0])
initial_x = inverse_transform(initial_params)
print("=" * 60)
print("REACTOR OPTIMIZATION PROBLEM")
print("=" * 60)
print("\nObjective: maximize (conversion - 0.01 × cost)")
print(" where cost = V + 0.001×(T-300)²")
print(f"\nConstraints:")
print(f" Volume: {V_bounds[0]} ≤ V ≤ {V_bounds[1]} m³")
print(f" Temperature: {T_bounds[0]} ≤ T ≤ {T_bounds[1]} K")
print(f"\nInitial guess: V = {initial_params[0]:.2f} m³, T = {initial_params[1]:.1f} K")
# Optimize using BFGS
solver = optx.BFGS(rtol=1e-5, atol=1e-5)
solution = optx.minimise(bounded_objective, solver, initial_x, max_steps=100, throw=False)
# Get optimal parameters
optimal_params = transform_params(solution.value)
V_opt, T_opt = float(optimal_params[0]), float(optimal_params[1])
X_opt = float(conversion_vs_VT(optimal_params))
cost_opt = V_opt + 0.001 * (T_opt - 300.0) ** 2
profit_opt = X_opt - 0.01 * cost_opt
print("\n" + "=" * 60)
print("OPTIMIZATION RESULTS")
print("=" * 60)
print(f"\nOptimal design:")
print(f" Volume V* = {V_opt:.3f} m³")
print(f" Temperature T* = {T_opt:.1f} K")
print(f"\nPerformance at optimum:")
print(f" Conversion = {X_opt*100:.1f}%")
print(f" Cost = {cost_opt:.2f}")
print(f" Profit = {profit_opt:.4f}")
# Interpret the solution
print("\n" + "-" * 60)
print("INTERPRETATION")
print("-" * 60)
# Check V bound
if abs(V_opt - V_bounds[1]) < 0.01:
print(f"\n• V = {V_opt:.1f} m³ is at the UPPER BOUND.")
print(" → More volume would improve profit, but we're space-limited.")
print(" → Consider expanding facility or adding a second reactor.")
elif abs(V_opt - V_bounds[0]) < 0.01:
print(f"\n• V = {V_opt:.1f} m³ is at the LOWER BOUND.")
print(" → Smaller reactor would be better, but not practical.")
else:
print(f"\n• V = {V_opt:.1f} m³ is an INTERIOR optimum.")
print(" → This balances conversion benefit vs capital cost.")
# Check T bound
if abs(T_opt - T_bounds[1]) < 0.5:
print(f"\n• T = {T_opt:.1f} K is at the UPPER BOUND.")
print(" → Higher temperature would help, but materials can't handle it.")
elif abs(T_opt - T_bounds[0]) < 0.5:
print(f"\n• T = {T_opt:.1f} K is at the LOWER BOUND.")
print(" → Lower temperature preferred to minimize energy cost.")
else:
print(f"\n• T = {T_opt:.1f} K is an INTERIOR optimum.")
print(" → This balances reaction rate benefit vs energy cost.")
print(" → Gradient ∂(profit)/∂T ≈ 0 at this point.")
============================================================
REACTOR OPTIMIZATION PROBLEM
============================================================
Objective: maximize (conversion - 0.01 × cost)
where cost = V + 0.001×(T-300)²
Constraints:
Volume: 0.1 ≤ V ≤ 5.0 m³
Temperature: 300.0 ≤ T ≤ 450.0 K
Initial guess: V = 0.50 m³, T = 320.0 K
============================================================
OPTIMIZATION RESULTS
============================================================
Optimal design:
Volume V* = 1.153 m³
Temperature T* = 331.5 K
Performance at optimum:
Conversion = 98.8%
Cost = 2.15
Profit = 0.9669
------------------------------------------------------------
INTERPRETATION
------------------------------------------------------------
• V = 1.2 m³ is an INTERIOR optimum.
→ This balances conversion benefit vs capital cost.
• T = 331.5 K is an INTERIOR optimum.
→ This balances reaction rate benefit vs energy cost.
→ Gradient ∂(profit)/∂T ≈ 0 at this point.
7. Validating AD Against Finite Differences#
Let’s verify that automatic differentiation gives the same answer as finite differences - but is exact!
def F_B_out(V: Array) -> Array:
"""Outlet F_B as function of reactor volume."""
params = CSTRParams(
V=V,
rate_fn=rate_function,
stoich=stoichiometry,
rate_params={"A": jnp.array(1e6), "Ea": jnp.array(50000.0)},
species_order=species_order,
dH_rxn=jnp.array([-50000.0]),
)
cstr = CSTR(params, thermo=thermo, mode="isothermal")
inlet = make_stream({"A": 10.0, "B": 0.0}, T=300.0, P=101325.0)
outlet, _ = cstr(inlet, T_spec=350.0)
return outlet["F_B"]
V = jnp.array(1.0)
ad_grad = float(jax.grad(F_B_out)(V))
print("Comparing AD gradient with finite differences:")
print(f" AD gradient: {ad_grad:.8f}")
print("\n Finite difference approximations:")
for eps in [1e-2, 1e-4, 1e-6, 1e-8]:
fd_grad = (float(F_B_out(V + eps)) - float(F_B_out(V - eps))) / (2 * eps)
rel_error = abs(fd_grad - ad_grad) / abs(ad_grad) * 100
print(f" ε = {eps:.0e}: grad = {fd_grad:.8f}, error = {rel_error:.4f}%")
print("\n✓ AD gradients match finite differences but are exact!")
Comparing AD gradient with finite differences:
AD gradient: 0.05170608
Finite difference approximations:
ε = 1e-02: grad = 0.05171119, error = 0.0099%
ε = 1e-04: grad = 0.05170608, error = 0.0000%
ε = 1e-06: grad = 0.05170608, error = 0.0000%
ε = 1e-08: grad = 0.05170602, error = 0.0001%
✓ AD gradients match finite differences but are exact!
Summary#
This notebook demonstrated:
Analysis |
Function |
Use Case |
|---|---|---|
|
First derivative |
Sensitivity, optimization |
|
Jacobian matrix |
Multi-output sensitivity |
|
Second derivatives |
Curvature, Newton methods |
Key Takeaways:
AD provides exact gradients (not approximations)
Computing all gradients is efficient (one backward pass)
Gradients enable optimization and uncertainty propagation
Higher-order derivatives reveal curvature and diminishing returns