numerical_integration/
lib.rs

1//!  # Numerical Integration
2//!
3//!  Algorithms and traits for numerical approximation in the Rust language
4//!
5//!  The purpose of this crate is to provide a selection of common algorithms
6//!  for approximating differential equations and to unify them under a common API
7//!  in order to make the process of using and testing different integration schemes
8//!  as easy as possible.
9//!
10//!  # How to use
11//!
12//!  The primary entry point into the API is through the `Integrator`,
13//!  `VelIntegrator`, and `AdaptiveIntegrator` traits. Each one has an `init()`
14//!  method and `step()`. They are all designed in such a way so that all values
15//!  and state of the integrator are stored *externally* and passed in through method
16//!  arguments exclusively.
17//!
18//!  Each trait represents a different class of algorithms requiring different
19//!  kinds of data. The `Integrator` trait takes in the current state and time,
20//!  the timestep, and a closure/function that can compute the derivative of the
21//!  state vector. `VelIntegrator` does the same, but also requires a function
22//!  that computes the velocity in an admittedly convoluted way. And
23//!  `AdaptiveIntegrator` takes in a minimum error value instead of a time-step.
24//!
25//!  In addition to these traits are traits that are like the above but adapted to
26//!  not include generics in the function signature so that it can be used as in
27//!  `dyn` types.
28//!
29//!  To use, you can either work with the traits generally and pass in a particular
30//!  implementor, or you can just use the various algorithms directly.
31//!
32//!  At the moment, this crate includes Velocity Verlet and methods in the
33//!  Runge-Kutta family (including Euler and RK4), but others (such as the linear
34//!  multistep methods) could be added if enough people want them.
35//!
36//!  # Current state of the project
37//!
38//!  This project is currently in hiatus for now, and it will probably remain as such
39//!  unless enough people express interest in it to be continued. If that *does* happen
40//!  though, expect the API to have breaking changes, as I am not quite
41//!  satisfied with the current design and certain features I wish to add require
42//!  a slight redesign.
43
44
45extern crate maths_traits;
46
47use maths_traits::algebra::module_like::*;
48use maths_traits::analysis::metric::*;
49use maths_traits::analysis::real::*;
50
51type Eval<'a,R,S> = &'a dyn Fn(R,S) -> S;
52
53pub trait Integrator {
54    fn init<R:Real, S:VectorSpace<R>, F:Fn(R, S) -> S>(&self, state: S, _dt:R, _force: F) -> Box<[S]> {Box::new([state])}
55    fn step<R:Real, S:VectorSpace<R>, F:Fn(R, S) -> S>(&self, time:R, state: &mut [S], dt:R, force: F) -> S;
56}
57
58pub trait Integrates<R:Real, S:VectorSpace<R>> {
59    fn init(&self, state: S, _dt:R, _force:Eval<R,S>) -> Box<[S]> {Box::new([state])}
60    fn step(&self, time:R, state: &mut [S], dt:R, force:Eval<R,S>) -> S;
61}
62
63impl<I:Integrator, R:Real, S:VectorSpace<R>> Integrates<R,S> for I {
64    fn init(&self, state: S, dt:R, force:Eval<R,S>) -> Box<[S]> {
65        Integrator::init(self, state, dt, &*force)
66    }
67    fn step(&self, time:R, state: &mut [S], dt:R, force:Eval<R,S>) -> S {
68        Integrator::step(self, time, state, dt, &*force)
69    }
70}
71
72pub trait VelIntegrator {
73    fn init_with_vel<
74        R:Real, S:VectorSpace<R>, V:Fn(R,S)->S, F:Fn(R,S)->S
75    > (&self, state: S, _dt:R, _vel:V, _force:F) -> Box<[S]> {Box::new([state])}
76
77    fn step_with_vel<
78        R:Real, S:VectorSpace<R>, V:Fn(R,S)->S, F:Fn(R,S)->S
79    >(&self, time:R, state: &mut [S], dt:R, velocity:V, force:F) -> S;
80}
81
82pub trait VelIntegrates<R:Real, S:VectorSpace<R>> {
83    fn init_with_vel(&self, state: S, _dt:R, _vel:Eval<R,S>, _force:Eval<R,S>) -> Box<[S]> {Box::new([state])}
84    fn step_with_vel(&self, time:R, state: &mut [S], dt:R, vel:Eval<R,S>, force:Eval<R,S>) -> S;
85}
86
87impl<I:VelIntegrator, R:Real, S:VectorSpace<R>> VelIntegrates<R,S> for I {
88    fn init_with_vel(&self, state: S, dt:R, vel:Eval<R,S>, force:Eval<R,S>) -> Box<[S]> {
89        VelIntegrator::init_with_vel(self, state, dt, &*vel, &*force)
90    }
91    fn step_with_vel(&self, time:R, state: &mut [S], dt:R, vel:Eval<R,S>, force:Eval<R,S>) -> S {
92        VelIntegrator::step_with_vel(self, time, state, dt, &*vel, &*force)
93    }
94}
95
96pub trait AdaptiveIntegrator {
97    fn adaptive_init<R:Real, S:VectorSpace<R>, M:Metric<S,R>, F:Fn(R, S) -> S>(&self, t0:R, state: S, _ds:R, _force:F, _d:M) -> Box<[(R,S)]>{
98        Box::new([(t0, state)])
99    }
100    fn adaptive_step<R:Real, S:VectorSpace<R>, M:Metric<S,R>, F:Fn(R, S) -> S>(&self, state: &mut [(R,S)], ds:R, force:F, d:M) -> (R,S);
101}
102
103pub use runge_kutta::*;
104pub mod runge_kutta;
105
106#[derive(Clone, Copy, PartialEq, Eq, Debug, Hash)]
107pub struct VelocityVerlet;
108
109impl VelIntegrator for VelocityVerlet {
110    fn init_with_vel<R:Real, S:VectorSpace<R>, V:Fn(R,S)->S, F:Fn(R,S)->S>(
111        &self, state: S, _dt:R, _vel:V, _force:F
112    ) -> Box<[S]> {
113        Box::new([state, S::zero()])
114    }
115
116    fn step_with_vel<R:Real, S:VectorSpace<R>, V:Fn(R,S)->S, F:Fn(R,S)->S>(
117        &self, time:R, state: &mut [S], dt:R, velocity:V, force:F
118    ) -> S {
119
120        let (s1, rest) = state.split_first_mut().unwrap();
121        let (a1, _) = rest.split_first_mut().unwrap();
122
123        let mid =
124            s1.clone() +
125            a1.clone()*dt.clone() +
126            velocity(time.clone(), a1.clone()) * (dt.clone()*dt.clone()*R::repr(0.5));
127
128        *s1 += a1.clone()*(dt.clone()*R::repr(0.5));
129        *a1 = force(time.clone() + dt.clone(), mid);
130        *s1 += a1.clone()*(dt.clone()*R::repr(0.5));
131
132        s1.clone()
133    }
134}
numerical_integration/lib.rs

numerical_integration/
lib.rs
