Skip to content

Commit

Permalink
Merge pull request #27 from martinjrobins/i23-state-ownership
Browse files Browse the repository at this point in the history
I23 state ownership
  • Loading branch information
martinjrobins authored Apr 5, 2024
2 parents cacc3b6 + ac5991e commit 0e4da83
Showing 1 changed file with 14 additions and 84 deletions.
98 changes: 14 additions & 84 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,14 @@
# Diffsol

<img src="https://github.com/martinjrobins/diffsol/actions/workflows/rust.yml/badge.svg" alt="CI build status badge">
<div align="center">
<a href="https://docs.rs/diffsol">
<img src="https://img.shields.io/crates/v/diffsol.svg?label=docs&color=blue&logo=rust" alt="docs.rs badge">
</a>
<a href="https://github.com/martinjrobins/diffsol/actions/workflows/rust.yml">
<img src="https://github.com/martinjrobins/diffsol/actions/workflows/rust.yml/badge.svg" alt="CI build status badge">
</a>
</div>

# DiffSol

Diffsol is a library for solving ordinary differential equations (ODEs) or
semi-explicit differential algebraic equations (DAEs) in Rust. You can use it
Expand All @@ -22,20 +30,9 @@ Users can specify the equations to solve in the following ODE form:
M \frac{dy}{dt} = f(t, y, p)
```

where `M` is a mass matrix, `y` is the state vector, `t` is the time, `p` is a
vector of parameters, and `f` is the right-hand side function. The mass matrix
`M` is optional (assumed to be the identity matrix if not provided).

The RHS function `f` can be specified as a function
that takes the state vector `y`, the parameter vector `p`, and the time `t` as
arguments and returns `f(t, y, p)`. The action of the jacobian `J` of `f` must also be specified as a
function that takes the state vector `y`, the parameter vector `p`, the time `t`
and a vector `v` and returns the product `Jv`.

The action of the mass matrix `M` can also be specified as a function that takes an input vector `v`,
the parameter vector `p` and the time `t` and returns the product mass matrix `Mv`.
Note that the only requirement for this mass matrix operator is that it must be linear.
It can be, for example, a singular matrix with zeros on the diagonal (i.e. defining a semi-explicit DAE).
where $M$ is a mass matrix, $y$ is the state vector, $t$ is the time, $p$ is a
vector of parameters, and $f$ is the right-hand side function. The mass matrix
$M$ is optional (assumed to be the identity matrix if not provided).

## Installation

Expand All @@ -54,71 +51,4 @@ cargo add diffeq

## Usage

This example solves the Robertson chemical kinetics problem, which is a stiff ODE with the equations:

```math
\begin{align*}
\frac{dy_1}{dt} &= -0.04 y_1 + 10^4 y_2 y_3 \\
\frac{dy_2}{dt} &= 0.04 y_1 - 10^4 y_2 y_3 - 3 \times 10^7 y_2^2 \\
\frac{dy_3}{dt} &= 3 \times 10^7 y_2^2
\end{align*}
```

with the initial conditions `y_1(0) = 1`, `y_2(0) = 0`, and `y_3(0) = 0`. We set
the tolerances to `1.0e-4` for the relative tolerance and `1.0e-8`, `1.0e-6`,
and `1.0e-6` for the absolute tolerances of `y_1`, `y_2`, and `y_3`,
respectively. We set the problem up with the following code:

```rust
type T = f64;
type V = nalgebra::DVector<T>;
let problem = OdeBuilder::new()
.p([0.04, 1.0e4, 3.0e7])
.rtol(1e-4)
.atol([1.0e-8, 1.0e-6, 1.0e-6])
.build_ode(
|x: &V, p: &V, _t: T, y: &mut V| {
y[0] = -p[0] * x[0] + p[1] * x[1] * x[2];
y[1] = p[0] * x[0] - p[1] * x[1] * x[2] - p[2] * x[1] * x[1];
y[2] = p[2] * x[1] * x[1];
},
|x: &V, p: &V, _t: T, v: &V, y: &mut V| {
y[0] = -p[0] * v[0] + p[1] * v[1] * x[2] + p[1] * x[1] * v[2];
y[1] = p[0] * v[0]
- p[1] * v[1] * x[2]
- p[1] * x[1] * v[2]
- 2.0 * p[2] * x[1] * v[1];
y[2] = 2.0 * p[2] * x[1] * v[1];
},
|_p: &V, _t: T| V::from_vec(vec![1.0, 0.0, 0.0]),
)
.unwrap();

let mut solver = Bdf::default();
```

We can then solve the problem up to time `t = 0.4` with the following code:

```rust
let t = 0.4;
let _y = solver.solve(&problem, t).unwrap();
```

Or if we want to explicitly step through time and have access to the solution
state, we can use the following code:

```rust
let mut state = OdeSolverState::new(&problem);
solver.set_problem(&mut state, &problem);
while state.t <= t {
solver.step(&mut state).unwrap();
}
let _y = solver.interpolate(&state, t).unwrap();
```

Note that `step` will advance the state to the next time step as chosen by the
solver, and `interpolate` will interpolate the solution back to the exact time
`t` that is requested.



For more documentation and examples, see the [API documentation](https://docs.rs/diffsol/latest/diffsol/).

0 comments on commit 0e4da83

Please sign in to comment.