Crate level document

src/lib.rs

@@ -1,4 +1,65 @@
-//! Configurable ODE/PDE solver
+//! Configurable ODE solver
+//!
+//! Design
+//! -------
+//! When we try to solve initial value problem (IVP) of an ordinal differential equation (ODE),
+//! we have to specify
+//!
+//! - the model space. Writing the ODE as a form $dx/dt = f(x)$,
+//!   the linear space where $x$ belong to, e.g. $\mathbb{R}^n$ or $\mathbb{C}^n$
+//!   is called model space, and represented by [ModelSpec] trait.
+//! - the equation, i.e. $f$ of $dx/dt = f(x)$,
+//!   e.g. Lorenz three variable equation, single or multiple pendulum and so on.
+//! - the scheme, i.e. how to solve given ODE,
+//!   e.g. explicit Euler, Runge-Kutta, symplectic Euler, and so on.
+//!
+//! Some equations requires some schemes.
+//! For example, Hamilton systems require symplectic schemes,
+//! or stiff equations require semi- or full-implicit schemes.
+//! We would like to implement these schemes without fixing ODE,
+//! but its abstraction depends on each scheme.
+//! Explicit schemes assumes the ODE is in a form
+//! $$
+//! \frac{dx}{dt} = f(x, t)
+//! $$
+//! and hope to abstract $f$, but symplectic schemes assumes the ODE must be defined with Hamiltonian $H$
+//! $$
+//! \frac{\partial p}{\partial t} = -\frac{\partial H}{\partial q},
+//! \frac{\partial q}{\partial t} = \frac{\partial H}{\partial p}.
+//! $$
+//! Some equation may be compatible several abstractions.
+//! Hamiltonian systems can be integrated with explicit schemes
+//! by ignoring phase-space volume contraction,
+//! or stiff systems can be integrated with explicit schemes with very small time steps.
+//!
+//! This crate introduces traits for each abstractions, e.g. [Explicit] or [SemiImplicit],
+//! which are implemented for each equations corresponds to ODE itself, e.g. [ode::Lorenz63].
+//! Schemes, e.g. [explicit::Euler], use this traits as type-bound.
+//!
+//! <!--
+//! KaTeX auto render
+//! Crate-level proc-macro is currently unstable feature.
+//! -->
+//! <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/katex@0.13.13/dist/katex.min.css" integrity="sha384-RZU/ijkSsFbcmivfdRBQDtwuwVqK7GMOw6IMvKyeWL2K5UAlyp6WonmB8m7Jd0Hn" crossorigin="anonymous">
+//! <script defer src="https://cdn.jsdelivr.net/npm/katex@0.13.13/dist/katex.min.js" integrity="sha384-pK1WpvzWVBQiP0/GjnvRxV4mOb0oxFuyRxJlk6vVw146n3egcN5C925NCP7a7BY8" crossorigin="anonymous"></script>
+//! <script defer src="https://cdn.jsdelivr.net/npm/katex@0.13.13/dist/contrib/auto-render.min.js" integrity="sha384-vZTG03m+2yp6N6BNi5iM4rW4oIwk5DfcNdFfxkk9ZWpDriOkXX8voJBFrAO7MpVl" crossorigin="anonymous"></script>
+//! <script>
+//!     document.addEventListener("DOMContentLoaded", function() {
+//!         renderMathInElement(document.body, {
+//!           // customised options
+//!           // • auto-render specific keys, e.g.:
+//!           delimiters: [
+//!               {left: '$$', right: '$$', display: true},
+//!               {left: '$', right: '$', display: false},
+//!               {left: '\\(', right: '\\)', display: false},
+//!               {left: '\\[', right: '\\]', display: true}
+//!           ],
+//!           // • rendering keys, e.g.:
+//!           throwOnError : false
+//!         });
+//!     });
+//! </script>
+//!
 
 pub mod adaptor;
 pub mod explicit;

src/traits.rs

@@ -1,41 +1,102 @@
-//! Fundamental traits for solving EoM
-
 use ndarray::*;
 use ndarray_linalg::*;
 use num_traits::Float;
 
-/// Model specification
+#[cfg(doc)]
+use crate::{explicit::*, ode::*, semi_implicit::*};
+
+#[cfg_attr(doc, katexit::katexit)]
+/// Model space, the linear space where the system state is represented.
+///
+/// For an ODE in a form $dx/dt = f(x)$, the linear space where $x$ belongs to
+/// is the model space.
+/// It is usually $\mathbb{R}^n$ or $\mathbb{C}^n$,
+/// but this crate allows it in multi-dimensional e.g. $\mathbb{C}^{N_x \times N_y}$
+/// to support spectral methods for PDE whose state space is Fourier coefficients.
+///
 pub trait ModelSpec: Clone {
     type Scalar: Scalar;
     type Dim: Dimension;
+    /// Number of scalars to describe the system state.
     fn model_size(&self) -> <Self::Dim as Dimension>::Pattern;
 }
 
-/// Interface for time-step
+/// Interface for set/get time step for integration
 pub trait TimeStep {
     type Time: Scalar + Float;
     fn get_dt(&self) -> Self::Time;
     fn set_dt(&mut self, dt: Self::Time);
 }
 
-/// Core implementation for explicit schemes
+#[cfg_attr(doc, katexit::katexit)]
+/// Abstraction for implementing explicit schemes
+///
+/// Consider equation of motion of an autonomous system
+/// described as an initial value problem of ODE:
+/// $$
+/// \frac{dx}{dt} = f(x),\space x(0) = x_0
+/// $$
+/// where $x = x(t)$ describes the system state specified by [ModelSpec] trait at a time $t$.
+/// This trait specifies $f$ itself, i.e. this trait will be implemented for structs corresponds to
+/// equations like [Lorenz63],
+/// and used while implementing integrate algorithms like [Euler],
+/// [Heun], and [RK4] algorithms.
+/// Since these algorithms are independent from the detail of $f$,
+/// this trait abstracts this part.
+///
 pub trait Explicit: ModelSpec {
-    /// calculate right hand side (rhs) of Explicit from current state
+    /// Evaluate $f(x)$ for a given state $x$
+    ///
+    /// This requires `&mut self` because evaluating $f$ may require
+    /// additional internal memory allocated in `Self`.
+    ///
     fn rhs<'a, S>(&mut self, x: &'a mut ArrayBase<S, Self::Dim>) -> &'a mut ArrayBase<S, Self::Dim>
     where
         S: DataMut<Elem = Self::Scalar>;
 }
 
-/// Core implementation for semi-implicit schemes
+#[cfg_attr(doc, katexit::katexit)]
+/// Abstraction for implementing semi-implicit schemes for stiff equations
+///
+/// Consider equation of motion of a stiff autonomous system
+/// described as an initial value problem of ODE:
+/// $$
+/// \frac{dx}{dt} = Ax + f(x),\space x(0) = x_0
+/// $$
+/// where $x = x(t)$ describes the system state specified by [ModelSpec] trait.
+/// We split the right hand side of the equation
+/// as the linear part $Ax$ to be stiff and the nonlinear part $f(x)$ not to be stiff.
+/// In addition, we assume $A$ is diagonalizable,
+/// and $x$ is selected to make $A$ diagonal.
+/// Similar to [Explicit], this trait abstracts the pair $(A, f)$ to implement
+/// semi-implicit schemes like [DiagRK4].
+///
+/// Stiff equations and semi-implicit schemes
+/// -------------------------------------------
+/// The stiffness causes numerical instabilities.
+/// For example, consider solving one-dimensional ODE $dx/dt = -\lambda x$
+/// with large $\lambda$ using explicit Euler scheme.
+/// Apparently, the solution is $x(t) = x(0)e^{-\lambda t}$,
+/// which converges to $0$ very quickly.
+/// However, to capture this process using explicit scheme like Euler scheme,
+/// we need as small $\Delta t$ as $\lambda^{-1}$.
+/// Such small $\Delta t$ is usually unacceptable,
+/// and implicit schemes are used for stiff equations,
+/// but full implicit schemes require solving fixed point problem like
+/// $1 + \lambda f(x) = 0$, which makes another instabilities.
+/// Semi-implicit schemes has been introduced to resolve this situation,
+/// i.e. use implicit scheme only for stiff linear part $Ax$
+/// and use explicit schemes on non-stiff part $f(x)$.
+///
 pub trait SemiImplicit: ModelSpec {
-    /// non-linear part of stiff equation
+    /// Non-stiff part $f(x)$
     fn nlin<'a, S>(
         &mut self,
         x: &'a mut ArrayBase<S, Self::Dim>,
     ) -> &'a mut ArrayBase<S, Self::Dim>
     where
         S: DataMut<Elem = Self::Scalar>;
-    /// diagonal elements of stiff linear part
+    /// Diagonal elements of stiff linear part of $A$
     fn diag(&self) -> Array<Self::Scalar, Self::Dim>;
 }