differential-equations 0.6.1

A Rust library for solving differential equations.
Documentation
# Quick Start: Defining State for ODEs

This guide provides a brief overview of how to define the state and differential equations for use with the Ordinary Differential Equation (ODE) solvers.

## The `ODE` Trait

The core of defining your problem lies in implementing the `ODE` trait. This trait defines the differential equation `dydt = f(t, y)` that the solver will use.

Key components of the `ODE` trait:

*   `diff(t, y, dydt)`: This method defines the system of differential equations.
    *   `t`: The current value of the independent variable (e.g., time).
    *   `y`: A reference to the current state vector.
    *   `dydt`: A mutable reference to the vector where the derivatives should be stored.
*   `event(t, y)`: An optional method to define conditions for interrupting the solver.
    *   It receives the current time `t` and state `y`.
    *   It should return a `ControlFlag`. By default, it returns `ControlFlag::Continue`.
    *   To terminate, return `ControlFlag::Terminate`, where `reason` can be any type implementing `Clone + Debug` (defaults to `String`).

## Example Implementation

Here's an example of implementing the `ODE` trait for a logistic growth model:

```rust
// Includes required elements and common methods.
// Less common methods are in the modules such as`ode::methods::...`
use differential_equations::prelude::*; 

// Define a struct to hold parameters for the ODE
struct LogisticGrowth {
    k: f64, // Growth rate
    m: f64, // Carrying capacity
}

// Implement the ODE trait for the LogisticGrowth struct
impl ODE for LogisticGrowth {
    // Define the differential equation: dy/dt = k * y * (1 - y / m)
    fn diff(&self, _t: f64, y: &f64, dydt: &mut f64) {
        *dydt = self.k * *y * (1.0 - *y / self.m);
    }

    // Optional: Define an event function to stop the solver
    fn event(&self, _t: f64, y: &f64) -> ControlFlag {
        if *y > 0.9 * self.m { // If population exceeds 90% of carrying capacity
            ControlFlag::Terminate("Reached 90% of carrying capacity".to_string())
        } else {
            ControlFlag::Continue // Otherwise, continue solving
        }
    }
}
```

## Solving an Initial Value Problem with `IVP`

Once you have implemented the `ODE` trait for your system, use `IVP::ode` to set up
and solve the initial value problem. The `IVP` builder acts as a solver controller,
orchestrating the solution process.

### Setting up the Problem

You create an ODE IVP by providing:
*   An instance of your `ODE` implementation (e.g., `LogisticGrowth`).
*   The initial time `t0`.
*   The final time `tf`.
*   The initial state vector `y0`.

### Choosing a Numerical Method

You'll also need to select a numerical method (solver). The library offers various fixed-step and adaptive-step solvers. For example, `DOP853` is a good general-purpose adaptive solver.

### Solving and Handling Output

The `solve` method on the `IVP` builder executes the numerical integration. You can configure how the solution is outputted, for instance, by requesting points at even intervals using the `even()` method.

The `solve` method returns a `Result<Solution, Status>`. The `Solution` struct contains the time points, state vectors, and solver statistics. The `Status` indicates how the solver finished (e.g., `Ok`, `Interrupted`).

### Example: Solving Logistic Growth

Here's how to solve the `LogisticGrowth` model defined earlier:

```rust
use differential_equations::prelude::*; 

fn main() {
    // Choose a numerical method and set tolerances
    let method = ExplicitRungeKutta::dop853()
        .rtol(1e-7) // Set relative tolerance
        .atol(1e-7);// Set absolute tolerance

    // Define initial conditions and time span
    let y0 = [1.0];        // Initial population
    let t0 = 0.0;          // Start time
    let tf = 10.0;         // End time

    // Create an instance of the ODE
    let ode_system = LogisticGrowth { k: 1.0, m: 10.0 };

    // Solve the problem, requesting output every 1.0 time unit
    match IVP::ode(&ode_system, t0, tf, y0)
        .even(1.0) // Set Solout to output every 1.0 time unit
        .method(method)
        .solve() // Solve the ODE
    {
        Ok(solution) => {
            // Check if the solver was interrupted by the event function
            if let Status::Interrupted(ref reason) = solution.status {
                println!("Solver stopped: {}", reason);
            }

            // Print the solution
            println!("Solution (t, y):");
            for (t, y_val) in solution.iter() {
                println!("({:.4}, {:.4})", t, y_val[0]);
            }

            // Print solver statistics
            println!("\nSolver Statistics:");
            println!("  Function evaluations: {}", solution.evals.function);
            println!("  Steps taken: {}", solution.steps.total());
            println!("  Accepted steps: {}", solution.steps.accepted);
            println!("  Rejected steps: {}", solution.steps.rejected);
        }
        Err(e) => eprintln!("An error occurred: {:?}", e),
    }
}
```

### Output

```sh
Solver stopped: Reached 90% of carrying capacity
Solution (t, y):
(0.0000, 1.0000)
(1.0000, 2.3197)
(2.0000, 4.5085)
(3.0000, 6.9057)
(4.0000, 8.5849)
(4.3944, 9.0000)

Solver Statistics:
  Function evaluations: 359
  Steps taken: 25
  Accepted steps: 21
  Rejected steps: 4
```

## Understanding Generics: `T` and `Y`

The `ODE` trait and related components like `IVP::ode` are defined with generics `<T, Y,...>`:

*   `T`: Represents the floating-point type used for calculations (e.g., `f64` or `f32`). This type applies to time and the components of the state vector.
*   `Y`: Represents the type of the state vector. This can be:
    *   A scalar `T` (e.g., `f64` or `f32`) for single-component systems.
    *   A `[T; N]` fixed-size array.
    *   A `Vec<T>` dynamically sized vector.
    *   An `nalgebra::SVector<T, N>` or dynamic nalgebra matrix/vector with the `nalgebra` feature.
    *   An `ndarray::Array<T, D>` with the `ndarray` feature.
    *   A `faer::Mat<T>` with the `faer` feature.
    *   A custom struct (deriving the `State` trait) with fields of type `T`.

The examples directory intentionally uses a mix of these forms: scalars for the
smallest examples, `Vec<f64>` for dynamic vector states, `ndarray::Array1`
for an event example, `faer::Mat` for a matrix-backed scalar integration example,
nalgebra vectors/matrices for linear algebra-heavy examples, and custom structs
with `#[derive(State)]` for named domain state.

Omitted `ODE` generics default to `T = f64` and `Y = f64`.

In the `LogisticGrowth` example, `impl ODE` uses the default scalar
state. For multi-component systems, use `[f64; N]`, `Vec<f64>`, a feature-backed
matrix/array type, or a custom `State` struct.

This generic setup provides flexibility, allowing you to define simple or complex systems of differential equations tailored to your specific numerical precision and state representation needs.