# Radau User Guide
This is the guide for the maintained Radau implementation in
`src/numerical/Radau/new`. The archived implementation is kept only as a
reference and is not the default production route.
## 1. Choose The API
There are two public entry points.
`UniversalODESolver::radau(...)` is the small convenience constructor. It uses
the maintained fixed-order Radau route with symbolic Lambdify, ExprLegacy, and
dense storage. It is useful when one facade must switch between several
solvers and no Radau-specific policy is needed.
For frontend, layout, continuation, output, telemetry, execution policy, or
AOT control, use the direct API:
```rust
use RustedSciThe::numerical::Radau::{
RadauConfig, RadauFrontend, RadauMatrixLayout, RadauProblem, RadauSolver,
RadauTelemetryMode,
};
use RustedSciThe::symbolic::symbolic_engine::Expr;
let problem = RadauProblem::new(
vec![Expr::parse_expression("-k*y")],
vec!["y".to_owned()],
"t",
)
.with_parameters(vec!["k".to_owned()])
.with_jacobian(vec![Expr::parse_expression("-k")]);
let config = RadauConfig {
frontend: RadauFrontend::ExprLegacy,
matrix_layout: RadauMatrixLayout::Dense,
telemetry: RadauTelemetryMode::Timings,
t_bound: 1.0,
..RadauConfig::default()
};
let mut solver = RadauSolver::prepare(problem, config)?;
let solution = solver.solve_with_parameters(&[1.0], &[2.0])?;
println!("t={} y={:?}", solution.t, solution.y);
# Ok::<(), Box<dyn std::error::Error>>(())
```
Preparation is fallible and separate from solving. A prepared solver owns the
symbolic callbacks and can be reused for a parameter series.
## 2. Select A Frontend
`RadauFrontend::ExprLegacy` keeps symbolic preparation in the legacy expression
representation. It is a good small-system baseline and a useful independent
correctness reference.
`RadauFrontend::AtomViewNative` converts the model once to AtomView and keeps
the residual/Jacobian preparation and callback path native to that
representation. It is usually the first route to try for large symbolic
systems, but there is no universal winner.
Frontend choice is independent from runtime execution:
- `Lambdify` prepares in-process closures and has no external compiler setup.
- `Aot` generates, compiles, links, and reuses a native artifact.
Both frontends use the same Radau numerical core. Compare them on the same
equations, initial state, tolerances, layout, and output policy.
## 3. Select A Matrix Layout
Choose the layout once during preparation:
```rust
let dense = RadauMatrixLayout::Dense;
let sparse = RadauMatrixLayout::Sparse;
let banded = RadauMatrixLayout::Banded { lower: 1, upper: 1 };
```
Use `Dense` for fully coupled Jacobians and small or medium systems. Use
`Banded` when the Jacobian has a stable narrow bandwidth. Use `Sparse` when a
stable sparse pattern is available and the system is large enough for sparse
storage to amortize its bookkeeping.
There is no silver bullet. A diffusion chain often favors Banded or Sparse;
small combustion or fully coupled kinetics can favor Dense. The layout must
match the actual Jacobian structure. Do not select a structured layout merely
because the residual is cheap.
The backend is not silently converted between Dense, Sparse, and Banded during
scalar callback work. A prepared layout owns its corresponding callback and
linear workspace.
## 4. Jacobian Contract
For symbolic Radau, provide an explicit Jacobian whenever possible. The vector
is row-major: for `n` equations and `n` variables it contains
`df[0]/dy[0] .. df[0]/dy[n-1]`, then the second row, and so on.
```rust
let problem = RadauProblem::new(
vec![Expr::parse_expression("-k*y")],
vec!["y".to_owned()],
"t",
)
.with_parameters(vec!["k".to_owned()])
.with_jacobian(vec![Expr::parse_expression("-k")]);
```
`RadauJacobianSource::Analytic` is the default and is the recommended policy
for stiff production workloads. `FiniteDifference` is useful as a diagnostic
or fallback policy, but it adds residual evaluations and should be compared
against an analytic-Jacobian run. AOT and structured layouts should be tested
with an explicit Jacobian first; this makes the generated contract and pattern
validation unambiguous.
### Native callback API
When the model is already implemented as Rust callbacks, use
`RadauNativeSolver` instead of rebuilding a symbolic expression graph. The
residual callback is mandatory. Supplying a Jacobian selects the analytic
path; passing `None` selects the typed forward-difference fallback. This
boundary currently uses Dense matrices and does not silently convert a native
callback into Sparse or Banded storage.
```rust
use nalgebra::{DMatrix, DVector};
use RustedSciThe::numerical::Radau::{RadauConfig, RadauNativeSolver};
let mut analytic = RadauNativeSolver::prepare(
RadauConfig { t_bound: 0.5, ..RadauConfig::default() },
|_t, y: &DVector<f64>| DVector::from_vec(vec![-y[0]]),
Some(|_t: f64, _y: &DVector<f64>| DMatrix::from_vec(1, 1, vec![-1.0])),
)?;
let solution = analytic.solve(&[1.0])?;
let mut finite_difference = RadauNativeSolver::prepare(
RadauConfig { t_bound: 0.5, ..RadauConfig::default() },
|_t, y: &DVector<f64>| DVector::from_vec(vec![-y[0]]),
Option::<fn(f64, &DVector<f64>) -> DMatrix<f64>>::None,
)?;
let fd_solution = finite_difference.solve(&[1.0])?;
# let _ = (solution, fd_solution);
# Ok::<(), Box<dyn std::error::Error>>(())
```
The public configuration reports `RadauExecution::NativeCallbacks` and either
`RadauJacobianSource::Analytic` or `FiniteDifference`. With telemetry enabled,
the latter reports its internal residual probes as
`finite_difference_probes`; those probes are not user-level residual calls.
The native callback boundary is currently Dense-only. Use the symbolic
`RadauProblem` API for Sparse or Banded layouts; it preserves the selected
structured backend without converting a native dense callback at runtime.
## 5. Tolerances And Step Controls
The main numerical controls are:
```rust
let config = RadauConfig {
rtol: 1e-6,
atol: 1e-9,
first_step: None,
max_step: f64::INFINITY,
max_steps: 1_000_000,
max_newton_iterations: 6,
max_retries: 12,
..RadauConfig::default()
};
```
Start with the defaults, then tighten `rtol` and `atol` together with a
problem-scaled validation criterion. Set `max_step` when a known physical
scale or rapidly changing forcing requires it. Set `first_step` only when the
default initial-step estimate is inappropriate. Increasing Newton or retry
budgets can hide a bad Jacobian or poor scaling, so inspect telemetry before
raising them.
## 6. Single Solve And Continuation
For a parameter-free model:
```rust
let solution = solver.solve(&[1.0])?;
```
For a model declared with `with_parameters`, bind values explicitly:
```rust
for k in [0.5, 1.0, 2.0, 4.0] {
let solution = solver.continue_with_parameters(&[1.0], &[k])?;
println!("k={k}: y(1)={:.9e}", solution.y[0]);
}
```
Continuation reuses the prepared symbolic model. It does not repeat symbolic
differentiation or frontend preparation merely because parameter values
changed. The first solve pays preparation cost; later solves are the relevant
warm path. For a fair benchmark report the first preparation separately from
the repeated continuation solves.
To reuse the prepared callbacks on a new interval:
```rust
solver.restart(0.0, 2.0)?;
let solution = solver.continue_with_parameters(&[1.0], &[2.0])?;
```
`restart` changes the interval, not the prepared model. Supply a new initial
state explicitly to the following solve.
## 7. Output And Post-Processing
Output retention is a deliberate memory policy:
```rust
use RustedSciThe::numerical::Radau::RadauOutputPolicy;
let final_only = RadauOutputPolicy::FinalOnly;
let dense_output = RadauOutputPolicy::Dense;
let selected_times = RadauOutputPolicy::Sampled(vec![0.25, 0.5, 0.75, 1.0]);
```
`FinalOnly` is the lowest-overhead choice. `Dense` retains continuous-output
segments and enables `solution.sample(t)` and `solution.sample_many(×)`.
`Sampled` evaluates only requested output times. Do not enable dense output in
a large production run unless post-processing really needs intermediate
states.
## 8. Telemetry
Telemetry is disabled by default and is designed not to charge the normal
path when disabled:
```rust
let config = RadauConfig {
telemetry: RadauTelemetryMode::Counters,
..RadauConfig::default()
};
```
Use `Counters` for calls, accepted/rejected steps, Newton iterations,
factorizations, and linear solves. Use `Timings` when investigating a staging
or performance question. The report separates preparation and runtime scopes;
parent and child timings are diagnostic and must not be added together.
Important fields include residual/Jacobian calls, factorization and solve
counts, callback timing, Newton timing, linear timing, output timing, parameter
rebinds, and AOT build/link/publication timing. Record the frontend, execution
family, layout, output policy, and parameter phase alongside the numbers.
## 9. AOT Lifecycle
AOT is selected explicitly and requires an AOT configuration:
```rust
use RustedSciThe::numerical::Radau::{RadauAotConfig, RadauExecution};
let config = RadauConfig {
execution: RadauExecution::Aot,
frontend: RadauFrontend::AtomViewNative,
matrix_layout: RadauMatrixLayout::Dense,
aot: Some(RadauAotConfig::build_if_missing_release(
"target/radau-artifacts",
)),
..RadauConfig::default()
};
```
The supported lifecycle helpers are:
- `build_if_missing_release(path)`: compile only when the artifact is absent;
- `require_prebuilt()`: fail fast if the artifact is not already published;
- `rebuild_always_release(path)`: force a fresh build for cold lifecycle tests.
Use `BuildIfMissing` for normal development, `RequirePrebuilt` for deployment
or process-isolated consumers, and `RebuildAlways` only for cold preparation
measurements. Keep AOT build/link time separate from warm solve time.
## 10. Parallel Policy
The default is safe sequential execution:
```rust
use RustedSciThe::numerical::Radau::RadauExecutionPolicy;
let sequential = RadauExecutionPolicy::Sequential;
let forced = RadauExecutionPolicy::Parallel { min_work: 256 };
let automatic = RadauExecutionPolicy::Auto { min_work: 256 };
```
Use `Parallel` only after measuring the workload. `Auto` is intended to avoid
parallel dispatch when worker overhead cannot amortize. Compare callback-only
and full-solve measurements separately; a callback win can disappear under
linear algebra or controller costs.
## 11. Error Handling And Examples
Preparation and solving return `Result<_, RadauError>`. Match error kinds in
application code rather than parsing display strings. Shape errors, invalid
configuration, unsupported routes, callback failures, Newton failures, step
budget exhaustion, and output-range errors are distinct failure classes.
Runnable examples in this repository:
- `examples/radau_lambdify_single.rs`: one parameterized Lambdify solve;
- `examples/radau_lambdify_continuation.rs`: prepared-model continuation;
- `examples/radau_public_api_guide.rs`: AtomView, counters, and dense output;
- `examples/radau_backends_guide.rs`: frontend/layout comparison;
- `examples/radau_aot_single.rs` and `examples/radau_aot_continuation.rs`:
explicit AOT lifecycle examples;
- `examples/universal_ode_example.rs`: the compact universal facade.
Validate correctness on the same workload before comparing performance. Keep
separate baselines for preparation, callback-only execution, full solve, and
continuation. A result that is faster on one workload or layout is not a
general claim about Radau.