Skip to content

Commit

Permalink
Removed redundant code and merged into one RungeKutta struct. Cleaned…
Browse files Browse the repository at this point in the history
… up the code a bit
  • Loading branch information
ysar committed Mar 30, 2024
1 parent 9838e3a commit ee57865
Show file tree
Hide file tree
Showing 6 changed files with 240 additions and 639 deletions.
5 changes: 3 additions & 2 deletions Cargo.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[package]
name = "lazyivy"
version = "0.2.0"
version = "0.3.0"
edition = "2021"
license = "BSD-2-Clause-Patent"
description = "Lazy Runge-Kutta integration for initial value problems."
Expand All @@ -9,7 +9,8 @@ repository = "https://github.com/ysar/lazyivy"
exclude = [
"examples/*",
]

keywords = ["runge-kutta", "ode", "integration", "lazy-evaluation"]
categories = ["mathematics", "science", "aerospace"]


# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
Expand Down
95 changes: 50 additions & 45 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,54 +4,54 @@
[![Build](https://github.com/ysar/lazyivy/actions/workflows/build.yml/badge.svg)](https://github.com/ysar/lazyivy/actions/workflows/build.yml)
[![Documentation](https://img.shields.io/docsrs/lazyivy/latest)](https://docs.rs/lazyivy/latest/lazyivy/)

lazyivy is a Rust crate that provides tools to solve initial value problems of the form
`dY/dt = F(t, y)` using Runge-Kutta methods.
lazyivy is a Rust crate that provides tools to solve initial value problems of
the form `dY/dt = F(t, y)` using Runge-Kutta methods.

Fixed step-size Runge-Kutta algorithms are implemented using the the structs,
- [`RungeKutta`](explicit_single::RungeKutta) to solve an ODE of one variable, and
- [`RungeKuttaSystem`](explicit_system::RungeKuttaSystem) to solve a system of ODEs.
The algorithms are implemented using the struct `RungeKutta`, that implements
`Iterator`. The following Runge-Kutta methods are implemented currently, and
more will be added in the near future.
- Euler 1
- Ralston 2
- Huen-Euler 2(1)
- Bogacki-Shampine 3(2)
- Fehlberg 4(5)
- Dormand-Prince 5(4)

These include the following RK methods - Euler, Ralston
Where `p` is the order of the method and `(p*)` is the order of the embedded
error estimator, if it is present.

Also implemented are embedded Runge-Kutta methods that consider an adaptive step-size based on
the local error estimate from a lower order method. These are implemented using the structs,
- [`RungeKuttaAdaptive`](explicit_single::RungeKuttaAdaptive) to solve an ODE of one variable
using an adaptive step size.
- [`RungeKuttaSystemAdaptive`](explicit_system::RungeKuttaSystemAdaptive) to solve a system of ODEs
with an adaptive step size.
## Lazy integration
`RungeKutta` implements the `Iterator` trait. Each `.next()` call advances the
iteration to the next Runge-Kutta *step* and returns a tuple `(t, y)`, where
`t` is the dependent variable and `y` is `Array1<f64>`.

These include the following RK methods of order `p(p*)` - Fehlberg 4(5), Hueneuler, Dormand-Prince 4(5).
Note that each Runge-Kutta *step* contains `s` number of internal *stages*.
Using lazyivy, there is no way at present to access the integration values for
these inner stages. The `.next()` call returns the final result for each step,
summed over all stages.

Where `p` is the order of the method and `p*` is the order of the error estimator step.

All integration structs implement the `Iterator` trait. Each `.next()` call advances the iteration
to the next Runge-Kutta *step* and returns a tuple `(t, y)`, where `t` is the dependent variable and
`y` can be either `f64`, as in the case of `RungeKutta` or `Array1<f64>`, in the case of
`RungeKuttaSystem`.

Note that each Runge-Kutta *step* contains `s` number of internal *stages*. Using lazyivy, there is
no way at present to access the integration values for these inner stages. The `next()` call returns
the final result for each step, summed over all stages.

The lazy implementation of Runge-Kutta means that you can consume the iterator in different ways.
For e.g., you can use `.last()` to keep only the final result, `.collect()` to gather the
state at all steps, `.map()` to chain the iterator with another, etc. You may also choose to use it
in a `for` loop and implement you own logic for modifying the step-size or customizing the stop
The lazy implementation of Runge-Kutta means that you can consume the iterator
in different ways. For e.g., you can use `.last()` to keep only the final
result, `.collect()` to gather the state at all steps, `.map()` to chain the
iterator with another, etc. You may also choose to use it in a `for` loop and
implement you own logic for modifying the step-size or customizing the stop
condition.

**API is still under development and future changes may not be backwards compatible.**
**API is unstable. It is active and under development.**

## Usage:

After adding lazyivy to `Cargo.toml`, create an initial value problem using
the various `new_*` methods. Here is an example
showing how to solve the [Brusselator](https://en.wikipedia.org/wiki/Brusselator).
```math

```math
\frac{d}{dt} \left[ \begin{array}{c}
y_1 \\ y_2 \end{array}\right] = \left[\begin{array}{c}1 - y_1^2 y_2 - 4 y_1 \\ 3y_1 - y_1^2 y_2 \end{array}\right]
y_1 \\ y_2 \end{array}\right] = \left[\begin{array}{c}1 - y_1^2 y_2 - 4 y_1
\\ 3y_1 - y_1^2 y_2 \end{array}\right]
```
```rust
use lazyivy::RungeKuttaSystemAdaptive;
use lazyivy::RungeKutta;
use ndarray::{Array, Array1};


Expand All @@ -68,32 +68,37 @@ fn main() {
let absolute_tol = Array::from_vec(vec![1.0e-4, 1.0e-4]);
let relative_tol = Array::from_vec(vec![1.0e-4, 1.0e-4]);

// Instantiate a integrator for an ODE system with adaptive step-size Runge-Kutta.
// Instantiate a integrator for an ODE system with adaptive step-size
// Runge-Kutta.

let mut integrator = RungeKuttaSystemAdaptive::new_fehlberg(
t0, // Initial condition - time
y0, // Initial condition - Brusselator variables in Array[y1, y2]
brusselator, // Evaluation function
|t, _| t > &20., // Predicate that determines stop condition
0.025, // Initial step size
relative_tol, // Relative tolerance for error estimation
absolute_tol, // Absolute tolerance for error estimation
let mut integrator = RungeKutta::new_fehlberg(
t0, // Initial condition - time t0
y0, // Initial condition - Initial condition [y1, y2] @ t0
brusselator, // Evaluation function
|t, _| t > &20., // Predicate that determines stop condition
0.025, // Initial step size
relative_tol, // Relative tolerance for error estimation
absolute_tol, // Absolute tolerance for error estimation
true, // Use adaptive step-size
);

// For adaptive algorithms, you can use this to improve the initial guess for the step size.
// For adaptive algorithms, you can use this to improve the initial guess
// for the step size.
integrator.set_step_size(&integrator.guess_initial_step());

// Perform the iterations and print each state.
for item in integrator {
println!("{:?}", item) // Prints the tuple (t, array[y1, y2]) at each iteration
println!("{:?}", item) // Prints (t, array[y1, y2]) for each step.
}
}
```
The result when plotted looks like this -
![Brusselator](https://raw.githubusercontent.com/ysar/lazyivy/main/examples/brusselator_adaptive.png)

To-do list:
## To-do list:
- [ ] Add more Runge-Kutta methods.
- [ ] Improve tests.
- [ ] Benchmark.
- [ ] Move allocations out of `next`.
- [x] Move allocations out of `next` and into a separate struct.
- [ ] Add more examples, e.g. Lorentz attractor.
- [ ]

0 comments on commit ee57865

Please sign in to comment.