Module peroxide::numerical::ode

source · [−]
Expand description

Solver for ordinary differential equations

Introduce ODE Trait & Structure

ODE Trait

  • ODE structures are divided by two kinds

    • ExplicitODE
    • ImplicitODE

  • ODE trait is given as

    extern crate peroxide;
use peroxide::fuga::{Real, State, BoundaryCondition, Environment};

pub trait ODE<E: Environment> {
    type Records;
    type Vector;
    type Param;
    type ODEMethod;

    fn mut_update(&mut self);
    fn integrate(&mut self) -> Self::Records;
    fn set_initial_condition<T: Real>(&mut self, init: State<T>) -> &mut Self;
    fn set_boundary_condition<T: Real>(
        &mut self,
        bound1: (State<T>, BoundaryCondition),
        bound2: (State<T>, BoundaryCondition),
    ) -> &mut Self;
    fn set_step_size(&mut self, dt: f64) -> &mut Self;
    fn set_method(&mut self, method: Self::ODEMethod) -> &mut Self;
    fn set_stop_condition(&mut self, f: fn(&Self) -> bool) -> &mut Self;
    fn set_times(&mut self, n: usize) -> &mut Self;
    fn check_enough(&self) -> bool;
    fn set_env(&mut self, env: E) -> &mut Self;
}
    • Records : The type to save results of ODE. Usually Matrix is used.
    • Vector : Vector can be below things.
      • Vec<f64> : Used for ExplicitODE
      • Vec<AD> : Used for ImplicitODE
    • Param : Also it can be f64 or AD
    • ODEMethod : Method for solving ODE
      • ExMethod : Explicit method
        • Euler : Euler first order
        • RK4 : Runge Kutta 4th order
      • ImMethod : Implicit method
        • BDF : Backward Euler 1st order (To be fixed)
        • GL4 : Gauss Legendre 4th order
    • Environment : External environment (CubicSpline, Vec, Matrix or Another external table)

State<T> structure

  • To use ODE trait, you should understand State<T> first.

    extern crate peroxide;
use peroxide::fuga::Real;

#[derive(Debug, Clone, Default)]
pub struct State<T: Real> {
    pub param: T,
    pub value: Vec<T>,
    pub deriv: Vec<T>,
}
    • T can be f64 or AD
      • f64 for ExplicitODE
      • AD for ImplicitODE
    • param is parameter for ODE. Usually it is represented by time.
    • value is value of each node.
    • deriv is value of derivative of each node.

For example,

$$ \frac{dy_n}{dt} = f(t, y_n) $$

  • $t$ is param
  • $y_n$ is value
  • $f(t,y_n)$ is deriv

Methods for State<T> are as follows.

  • to_f64(&self) -> State<f64>
  • to_ad(&self) -> State<AD>
  • new(T, Vec<T>, Vec<T>) -> Self

Environment

  • Environment needs just Default
  • To use custom Environment, just type follows : impl Environment for CustomType {}
  • If you don’t want to use Environment, then use NoEnv
  • Implemented Data Types
    • Vec<f64>
    • Polynomial
    • Matrix
    • CubicSpline
    • NoEnv
#[macro_use]
extern crate peroxide;
use peroxide::fuga::*;

fn main() {
    let x = seq(0, 10, 1);
    x.print();
    let y = x.iter().enumerate().map(|(i, &t)| t.powi(5-i as i32)).collect::<Vec<f64>>();

    let c = CubicSpline::from_nodes(&x, &y);

    let init_state = State::<f64>::new(0f64, c!(1), c!(0));

    let mut ode_solver = ExplicitODE::new(test_fn);

    ode_solver
        .set_method(ExMethod::RK4)
        .set_initial_condition(init_state)
        .set_step_size(0.01)
        .set_times(1000)
        .set_env(c);

    let result = ode_solver.integrate();
    result.print();
}

fn test_fn(st: &mut State<f64>, env: &CubicSpline) {
    let x = st.param;
    let dy = &mut st.deriv;
    dy[0] = env.eval(x);
}

ExplicitODE struct

ExplicitODE is given as follow :

extern crate peroxide;
use std::collections::HashMap;
use peroxide::fuga::{State, ExMethod, BoundaryCondition, ODEOptions, Environment};

#[derive(Clone)]
pub struct ExplicitODE<E: Environment> {
    state: State<f64>,
    func: fn(&mut State<f64>, &E),
    step_size: f64,
    method: ExMethod,
    init_cond: State<f64>,
    bound_cond1: (State<f64>, BoundaryCondition),
    bound_cond2: (State<f64>, BoundaryCondition),
    stop_cond: fn(&Self) -> bool,
    times: usize,
    to_use: HashMap<ODEOptions, bool>,
    env: E,
}
  • state : Current param, value, derivative
  • func : Function to update state
  • init_cond : Initial condition
  • bound_cond1 : If boundary problem, then first boundary condition
  • bound_cond2 : second boundary condition
  • stop_cond : Stop condition (stop before times)
  • times : How many times do you want to update?
  • to_use : Just check whether information is enough
  • env : Environment

ImplicitODE struct

ImplicitODE is struct to solve stiff ODE with implicit methods.

extern crate peroxide;
use std::collections::HashMap;
use peroxide::fuga::{State, ImMethod, BoundaryCondition, ODEOptions, Environment, AD};

#[derive(Clone)]
pub struct ImplicitODE<E: Environment> {
    state: State<AD>,
    func: fn(&mut State<AD>, &E),
    step_size: f64,
    method: ImMethod,
    init_cond: State<AD>,
    bound_cond1: (State<AD>, BoundaryCondition),
    bound_cond2: (State<AD>, BoundaryCondition),
    stop_cond: fn(&Self) -> bool,
    times: usize,
    to_use: HashMap<ODEOptions, bool>,
    env: E,
}

Example

Lorenz Butterfly

extern crate peroxide;
use peroxide::fuga::*;

fn main() {
    // =========================================
    //  Declare ODE (Euler)
    // =========================================
    let mut ex_test = ExplicitODE::new(f);

    let init_state: State<f64> = State::new(
        0.0,
        vec![10.0, 1.0, 1.0],
        vec![0.0, 0.0, 0.0],
    );

    ex_test
        .set_initial_condition(init_state)
        .set_method(ExMethod::Euler)
        .set_step_size(0.01f64)
        .set_times(10000);

    // =========================================
    //  Declare ODE (RK4)
    // =========================================
    let mut ex_test2 = ex_test.clone();
    ex_test2.set_method(ExMethod::RK4);

    // =========================================
    //  Declare ODE (GL4)
    // =========================================
    let mut im_test = ImplicitODE::new(g);

    let init_state: State<AD> = State::new(
        AD0(0.0),
        vec![AD0(10f64), AD0(1f64), AD0(1f64)],
        vec![AD0(0.0), AD0(0.0), AD0(0.0)],
    );

    im_test
        .set_initial_condition(init_state)
        .set_method(ImMethod::GL4)
        .set_step_size(0.01f64)
        .set_times(10000);

    // =========================================
    //  Save results
    // =========================================
    let results = ex_test.integrate();
    let results2 = ex_test2.integrate();
    let results3 = im_test.integrate();

    // Extract data
    let mut df = DataFrame::new(vec![]);
    df.push("x_euler", Series::new(results.col(1)));
    df.push("z_euler", Series::new(results.col(3)));
    df.push("x_rk4", Series::new(results2.col(1)));
    df.push("z_rk4", Series::new(results2.col(3)));
    df.push("x_gl4", Series::new(results3.col(1)));
    df.push("z_gl4", Series::new(results3.col(3)));

    // Write netcdf file (`nc` feature required)
    df.write_nc("example_data/lorenz.nc")
        .expect("Can't write lorenz.nc");
}

fn f(st: &mut State<f64>, _: &NoEnv) {
    let x = &st.value;
    let dx = &mut st.deriv;
    dx[0] = 10f64 * (x[1] - x[0]);
    dx[1] = 28f64 * x[0] - x[1] - x[0] * x[2];
    dx[2] = -8f64/3f64 * x[2] + x[0] * x[1];
}

fn g(st: &mut State<AD>, _: &NoEnv) {
    let x = &st.value;
    let dx = &mut st.deriv;
    dx[0] = 10f64 * (x[1] - x[0]);
    dx[1] = 28f64 * x[0] - x[1] - x[0] * x[2];
    dx[2] = -8f64/3f64 * x[2] + x[0] * x[1];
}

If plotting pickle data with python, then

Lorenz with Euler

Lorenz with RK4

Lorenz with GL4

Simple 1D Runge-Kutta

$$\begin{gathered} \frac{dy}{dx} = \frac{5x^2 - y}{e^{x+y}} \\ y(0) = 1 \end{gathered}$$

#[macro_use]
extern crate peroxide;
use peroxide::fuga::*;

fn main() {
    let init_state = State::<f64>::new(0f64, c!(1), c!(0));

    let mut ode_solver = ExplicitODE::new(test_fn);

    ode_solver
        .set_method(ExMethod::RK4)
        .set_initial_condition(init_state)
        .set_step_size(0.01)
        .set_times(1000);

    let result = ode_solver.integrate();

    // Plot or Extract..
}

fn test_fn(st: &mut State<f64>, _: &NoEnv) {
    let x = st.param;
    let y = &st.value;
    let dy = &mut st.deriv;
    dy[0] = (5f64*x.powi(2) - y[0]) / (x + y[0]).exp();
}

Structs

ExplicitODE
ImplicitODE
NoEnv
State
State for ODE

Enums

BoundaryCondition
Kinds of Boundary Conditions
ExMethod
Explicit ODE Methods
ImMethod
ODEOptions
Options for ODE

Traits

Environment
ODE
ODE solver