# BVP Solver User Guide
This guide provides detailed instructions for using the Boundary Value Problem (BVP) solver in RustedSciThe with the `faer` linear algebra backend.
## Table of Contents
1. [Basic Setup](#basic-setup)
2. [Defining the ODE System](#defining-the-ode-system)
3. [Setting Boundary Conditions](#setting-boundary-conditions)
4. [Initial Guess](#initial-guess)
5. [Solver Parameters](#solver-parameters)
6. [Complete Examples](#complete-examples)
7. [Advanced Features](#advanced-features)
## Basic Setup
First, import the necessary types and functions:
```rust
use crate::numerical::BVP_sci::BVP_sci_faer::solve_bvp;
use faer::col::Col;
use faer::mat::Mat;
type faer_col = Col<f64>;
type faer_dense_mat = Mat<f64>;
```
## Defining the ODE System
The ODE system is defined as a function that takes three parameters:
- `x`: mesh points (faer_col)
- `y`: solution values at mesh points (faer_dense_mat, shape: n×m)
- `p`: parameters (faer_col)
Returns: derivatives (faer_dense_mat, shape: n×m)
### Example: Second-order ODE y'' = -y
```rust
```rust
// y'' = y^3 - sin(x)
let fun = |x: &faer_col, y: &faer_dense_mat, _p: &faer_col| {
let mut f = faer_dense_mat::zeros(2, y.ncols());
for j in 0..y.ncols() {
let y1 = *y.get(0, j);
*f.get_mut(0, j) = *y.get(1, j);
*f.get_mut(1, j) = y1 * y1 * y1 - x[j].sin();
}
f
};
```
#### Parametric System
```rust
// y'' + p[0]*y' + p[1]*y = 0
let fun = |_x: &faer_col, y: &faer_dense_mat, p: &faer_col| {
let mut f = faer_dense_mat::zeros(2, y.ncols());
for j in 0..y.ncols() {
*f.get_mut(0, j) = *y.get(1, j);
*f.get_mut(1, j) = -p[1] * *y.get(0, j) - p[0] * *y.get(1, j);
}
f
};
```
## Setting Boundary Conditions
Boundary conditions are defined as a function that takes:
- `ya`: solution values at left boundary (faer_col)
- `yb`: solution values at right boundary (faer_col)
- `p`: parameters (faer_col)
Returns: boundary condition residuals (faer_col, length: n + k)
### Common Boundary Condition Types
#### Dirichlet Conditions
```rust
// y(0) = 1, y(1) = 0
```rust
// y(0) = y(1), y'(0) = y'(1)
let bc = |ya: &faer_col, yb: &faer_col, _p: &faer_col| {
faer_col::from_fn(2, |i| match i {
0 => ya[0] - yb[0], // y(0) = y(1)
1 => ya[1] - yb[1], // y'(0) = y'(1)
_ => 0.0,
})
};
```
#### Parameter Estimation
```rust
// y(0) = 0, y(1) = 1, ∫y dx = target_integral
let bc = |ya: &faer_col, yb: &faer_col, p: &faer_col| {
faer_col::from_fn(3, |i| match i {
0 => ya[0], // y(0) = 0
1 => yb[0] - 1.0, // y(1) = 1
2 => p[0] - target_integral, // parameter constraint
_ => 0.0,
})
};
```
## Initial Guess
The initial guess is a matrix of shape (n, m) where:
- n = number of equations
- m = number of mesh points
### Creating Initial Guess
#### Constant Values
```rust
let mut y = faer_dense_mat::zeros(2, 5);
for j in 0..5 {
*y.get_mut(0, j) = 1.0; // y1 = 1 everywhere
*y.get_mut(1, j) = 0.0; // y2 = 0 everywhere
}
```
#### Linear Interpolation
```rust
for j in 0..5 {
*y.get_mut(0, j) = x[j]; // linear from 0 to 1
*y.get_mut(1, j) = 1.0; // constant derivative
}
```
#### Analytical Approximation
```rust
// For y'' = -y with y(0)=0, y(π)=0, use sin(x) as initial guess
for j in 0..5 {
*y.get_mut(0, j) = x[j].sin();
*y.get_mut(1, j) = x[j].cos();
}
```
## Solver Parameters
### Basic Parameters
```rust
let result = solve_bvp(
&fun, // ODE function
&bc, // Boundary conditions
x, // Initial mesh
y, // Initial guess
None, // ODE Jacobian (optional)
None, // BC Jacobian (optional)
None, // Parameters (optional)
None, // Parameter bounds (optional)
1e-6, // BVP tolerance
1000, // Max iterations
2, // Max mesh refinements
None, // BC tolerance (optional)
);
```
### Parameter Descriptions
- **BVP tolerance** (`1e-6`): Controls accuracy of collocation residuals
- **Max iterations** (`1000`): Maximum Newton iterations per mesh
- **Max refinements** (`2`): Maximum number of adaptive mesh refinements
- **BC tolerance** (`None`): Boundary condition tolerance (defaults to BVP tolerance)
### Tolerance Guidelines
| Simple linear | 1e-8 | 1e-10 |
| Nonlinear smooth | 1e-6 | 1e-8 |
| Stiff/oscillatory | 1e-4 | 1e-6 |
| Parameter estimation | 1e-5 | 1e-7 |
## Complete Examples
### Example 1: Simple Linear BVP
```rust
use crate::numerical::BVP_sci::BVP_sci_faer::solve_bvp;
use faer::col::Col;
use faer::mat::Mat;
type faer_col = Col<f64>;
type faer_dense_mat = Mat<f64>;
fn solve_linear_bvp() {
// y'' = 0, y(0) = 0, y(1) = 1
// Exact solution: y = x
let fun = |_x: &faer_col, y: &faer_dense_mat, _p: &faer_col| {
let mut f = faer_dense_mat::zeros(2, y.ncols());
for j in 0..y.ncols() {
*f.get_mut(0, j) = *y.get(1, j); // y' = y2
*f.get_mut(1, j) = 0.0; // y2' = 0
}
f
};
let bc = |ya: &faer_col, yb: &faer_col, _p: &faer_col| {
faer_col::from_fn(2, |i| match i {
0 => ya[0], // y(0) = 0
1 => yb[0] - 1.0, // y(1) = 1
_ => 0.0,
})
};
// Initial mesh
let x = faer_col::from_fn(3, |i| i as f64 * 0.5);
// Initial guess
let mut y = faer_dense_mat::zeros(2, 3);
for j in 0..3 {
*y.get_mut(0, j) = x[j]; // linear guess
*y.get_mut(1, j) = 1.0; // constant derivative
}
let result = solve_bvp(
&fun, &bc, x, y, None, None, None, None,
1e-8, 1000, 0, None
);
match result {
Ok(res) => {
if res.success {
println!("Solution found!");
for i in 0..res.x.nrows() {
println!("x = {:.3}, y = {:.6}", res.x[i], res.y.get(0, i));
}
}
}
Err(e) => println!("Error: {}", e),
}
}
```
### Example 2: Nonlinear BVP with Parameters
```rust
fn solve_parametric_bvp() {
// y'' + p*y = 0, y(0) = 0, y(π) = 0, max(y) = 1
let fun = |_x: &faer_col, y: &faer_dense_mat, p: &faer_col| {
let mut f = faer_dense_mat::zeros(2, y.ncols());
for j in 0..y.ncols() {
*f.get_mut(0, j) = *y.get(1, j);
*f.get_mut(1, j) = -p[0] * *y.get(0, j);
}
f
};
let bc = |ya: &faer_col, yb: &faer_col, p: &faer_col| {
faer_col::from_fn(3, |i| match i {
0 => ya[0], // y(0) = 0
1 => yb[0], // y(π) = 0
2 => p[0] - 1.0, // p = 1 (parameter constraint)
_ => 0.0,
})
};
let x = faer_col::from_fn(5, |i| i as f64 * std::f64::consts::PI / 4.0);
let mut y = faer_dense_mat::zeros(2, 5);
for j in 0..5 {
*y.get_mut(0, j) = x[j].sin();
*y.get_mut(1, j) = x[j].cos();
}
let p = faer_col::from_fn(1, |_| 1.0); // Initial parameter guess
let result = solve_bvp(
&fun, &bc, x, y, None, None, Some(p), None,
1e-6, 1000, 2, None
);
match result {
Ok(res) => {
if res.success {
println!("Parameter found: p = {:.6}", res.p.as_ref().unwrap()[0]);
}
}
Err(e) => println!("Error: {}", e),
}
}
```
## Advanced Features
### Singular Term `S`
The solver now supports a SciPy-style singular-term matrix `S` on the
builder-based numerical route:
```rust
let options = NumericalBvpSolveOptions::new(mesh, guess, tol, max_nodes)
.with_singular_term(Some(singular_term));
```
This enables problems written in the form:
```text
dy/dx = f(x, y, p) + S * y / (x - a)
```
where `a` is the left endpoint of the mesh.
Practical notes:
- `S` must be square and match the state dimension.
- the current plumbing is wired into the numerical solve path first;
- if you are writing a closure-first workflow, `with_singular_term(...)`
belongs with the other solver options, not with symbolic setup.
### Custom Jacobians
For better performance and convergence, provide analytical Jacobians:
```rust
let fun_jac = |x: &faer_col, y: &faer_dense_mat, p: &faer_col| {
let m = x.nrows();
let mut df_dy = Vec::new();
let mut df_dp = Vec::new();
for j in 0..m {
// Jacobian w.r.t. y at each mesh point
let mut jac_y = faer::sparse::SparseColMat::try_new_from_triplets(
2, 2,
&[
faer::sparse::Triplet::new(0, 1, 1.0), // ∂f1/∂y2 = 1
faer::sparse::Triplet::new(1, 0, -p[0]), // ∂f2/∂y1 = -p
]
).unwrap();
df_dy.push(jac_y);
// Jacobian w.r.t. parameters
let mut jac_p = faer::sparse::SparseColMat::try_new_from_triplets(
2, 1,
&[faer::sparse::Triplet::new(1, 0, -*y.get(0, j))] // ∂f2/∂p = -y1
).unwrap();
df_dp.push(jac_p);
}
(df_dy, Some(df_dp))
};
3, 2,
&[faer::sparse::Triplet::new(0, 0, 1.0)] // ∂bc1/∂ya1 = 1
).unwrap();
let dbc_dyb = faer::sparse::SparseColMat::try_new_from_triplets(
3, 2,
&[faer::sparse::Triplet::new(1, 0, 1.0)] // ∂bc2/∂yb1 = 1
).unwrap();
let dbc_dp = faer::sparse::SparseColMat::try_new_from_triplets(
3, 1,
&[faer::sparse::Triplet::new(2, 0, 1.0)] // ∂bc3/∂p = 1
).unwrap();
(dbc_dya, dbc_dyb, Some(dbc_dp))
};
// Use with solver
let result = solve_bvp(
&fun, &bc, x, y,
Some(&fun_jac), // ODE Jacobian
Some(&bc_jac), // BC Jacobian
Some(p), None, 1e-6, 1000, 2, None
);
```
### Mesh Refinement Strategy
The solver automatically refines the mesh when:
- Collocation residuals are too large
- Solution changes rapidly between mesh points
- Maximum refinements not exceeded
Control refinement with the `max_refinements` parameter:
- `0`: No refinement (fixed mesh)
- `1-3`: Conservative refinement
- `4-10`: Aggressive refinement (may be slow)
### Error Handling
```rust
match solve_bvp(&fun, &bc, x, y, None, None, None, None, 1e-6, 1000, 2, None) {
Ok(result) => {
match result.status {
0 => println!("Success: Solution converged"),
1 => println!("Warning: Max refinements reached"),
2 => println!("Error: Singular Jacobian"),
3 => println!("Error: BC tolerance not met"),
_ => println!("Error: Max iterations reached"),
}
}
Err(e) => println!("Solver error: {}", e),
}
```
This guide covers the essential aspects of using the BVP solver. For more complex problems, consider starting with simpler initial guesses and gradually increasing accuracy requirements.
## Text task documents and backend selection
The task-document route is intended for reproducible examples and batch runs;
the typed native API remains the preferred application-facing interface. A
minimal BVP_sci task can make the solver and numerical route explicit:
```text
# Comments are ignored by the task parser.
solver: BVP_sci
method: Dense
frontend: ExprLegacy
execution: Lambdify
tolerance: 1e-6
max_iterations: 100
output_policy: none
```
`method` selects the matrix layout and `frontend` selects the symbolic route.
They are separate dimensions: `ExprLegacy` and `AtomViewNative` describe how
the residual and Jacobian are represented and prepared, while `Lambdify` and
`AOT` describe how the prepared callbacks are executed. The task parser also
accepts the BVP_Damp route as `solver: BVP`; do not confuse these two solver
identifiers when comparing numerical results.
There is no universal fastest matrix layout. Dense is a natural choice for
small fully coupled systems, while Sparse or Banded can be better when the
global collocation Jacobian has exploitable structure. Compare matched
workloads rather than comparing different equations. Preparation, callback,
linear solve, and full-solve timings must be reported separately: a faster
symbolic frontend does not imply a faster dense factorization.
For AOT task documents, use `aot_build_policy` with `BuildIfMissing`,
`RequirePrebuilt`, or `RebuildAlways`, together with `output_dir`, `compiler`,
and an optional publication/handoff directory. `RequirePrebuilt` is a strict
consumer mode and reports missing, stale, provenance, compiler, link, and
timeout failures as typed errors. The lifecycle record distinguishes cache
lookup, materialization, build, link, publication, and runtime readiness from
the numerical solve.
Continuation is expressed with a `continuation` block, for example:
```text
continuation:
parameter: k
values: 1.0, 1.5, 2.0, 2.5
mode: prepared
restart_each: false
```
`fresh` reparses and prepares every segment, `warm` reuses the prepared model
while rebinding parameters, and `prepared` requests the strongest reuse policy.
An explicit restart may provide a new `y0`, `t0`, and `t_end` for a segment.
Keep these policies distinct in performance reports and check memory retention
for long series. `output_policy` may be `none`, `terminal`, `plotters`, or
`gnuplot`; comments beginning with `#` are ignored. The shared task-document
grammar and batch-runner conventions are described in
[`TASK_DOCS_GUIDE_EN.md`](../../command_interpreter/TASK_DOCS_GUIDE_EN.md).