Expand description
Solver for ordinary differential equations
Introduce ODE
Trait & Structure
ODE
Trait
-
ODE
structures are divided by two kindsExplicitODE
ImplicitODE
-
ODE
trait is given asextern 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. UsuallyMatrix
is used.Vector
: Vector can be below things.Vec<f64>
: Used forExplicitODE
Vec<AD>
: Used forImplicitODE
Param
: Also it can bef64
orAD
ODEMethod
: Method for solving ODEExMethod
: Explicit methodEuler
: Euler first orderRK4
: Runge Kutta 4th order
ImMethod
: Implicit methodBDF
: 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 understandState<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 bef64
orAD
f64
forExplicitODE
AD
forImplicitODE
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 justDefault
- To use custom
Environment
, just type follows :impl Environment for CustomType {}
- If you don’t want to use
Environment
, then useNoEnv
- 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, derivativefunc
: Function to updatestate
init_cond
: Initial conditionbound_cond1
: If boundary problem, then first boundary conditionbound_cond2
: second boundary conditionstop_cond
: Stop condition (stop beforetimes
)times
: How many times do you want to update?to_use
: Just check whether information is enoughenv
: 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
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
Enums
Kinds of Boundary Conditions
Explicit ODE Methods
Options for ODE
Traits
ODE solver