150 lines
4.7 KiB
Rust
150 lines
4.7 KiB
Rust
#![allow(dead_code)]
|
|
|
|
pub mod callback;
|
|
pub mod controller;
|
|
pub mod integrator;
|
|
pub mod ode;
|
|
pub mod problem;
|
|
|
|
pub mod prelude {
|
|
pub use super::callback::{stop, Callback};
|
|
pub use super::controller::PIController;
|
|
pub use super::integrator::dormand_prince::DormandPrince45;
|
|
pub use super::ode::ODE;
|
|
pub use super::problem::{Problem, Solution};
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
use crate::prelude::*;
|
|
use approx::assert_relative_eq;
|
|
use nalgebra::{Vector1, Vector2, Vector6};
|
|
use std::f64::consts::PI;
|
|
|
|
#[test]
|
|
fn test_readme() {
|
|
// Define the system (parameters, derivative, and initial state)
|
|
type Params = (f64, f64); // Gravity and Length of Pendulum
|
|
let params = (9.81, 1.0);
|
|
|
|
fn derivative(_t: f64, y: Vector2<f64>, p: &Params) -> Vector2<f64> {
|
|
let &(g, l) = p;
|
|
let theta = y[0];
|
|
let d_theta = y[1];
|
|
Vector2::new(d_theta, -(g / l) * theta.sin())
|
|
}
|
|
|
|
let y0 = Vector2::new(0.0, PI / 2.0);
|
|
|
|
// Set up the problem (ODE, Integrator, Controller, and Callbacks)
|
|
let ode = ODE::new(&derivative, 0.0, 6.3, y0, params);
|
|
let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6);
|
|
let controller = PIController::default();
|
|
|
|
let value_too_high = Callback {
|
|
event: &|t: f64, _y: Vector2<f64>, _p: &Params| 5.0 - t,
|
|
effect: &stop,
|
|
};
|
|
|
|
// Solve the problem
|
|
let mut problem = Problem::new(ode, dp45, controller).with_callback(value_too_high);
|
|
let solution = problem.solve();
|
|
|
|
// Can interpolate solutions to whatever you want
|
|
let _interpolated_answer = solution.interpolate(4.4);
|
|
}
|
|
|
|
#[test]
|
|
fn test_correctness() {
|
|
// Define the system (parameters, derivative, and initial state)
|
|
type Params = ();
|
|
let params = ();
|
|
|
|
fn derivative(_t: f64, y: Vector1<f64>, _p: &Params) -> Vector1<f64> {
|
|
Vector1::new(5.0 * y[0] - 3.0)
|
|
}
|
|
|
|
let y0 = Vector1::new(1.0);
|
|
|
|
// Set up the problem (ODE, Integrator, Controller, and Callbacks)
|
|
let ode = ODE::new(&derivative, 2.0, 3.0, y0, params);
|
|
let dp45 = DormandPrince45::new();
|
|
let controller = PIController::default();
|
|
|
|
// Solve the problem
|
|
let mut problem = Problem::new(ode, dp45, controller);
|
|
let solution = problem.solve();
|
|
for (time, state) in solution.times.iter().zip(solution.states.iter()) {
|
|
let exact = 0.4 * (5.0 * (time - 2.0)).exp() + 0.6;
|
|
assert_relative_eq!(state[0], exact, max_relative = 1e-7);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn test_orbit() {
|
|
// Calculate one period
|
|
let a = 6.7781363e6_f64;
|
|
let mu = 3.98600441500000e14;
|
|
let period = 2.0 * PI * (a.powi(3) / mu).sqrt();
|
|
|
|
// Set up the system
|
|
type Params = (f64,);
|
|
let params = (mu,);
|
|
fn derivative(_t: f64, state: Vector6<f64>, p: &Params) -> Vector6<f64> {
|
|
let acc = -(p.0 * state.fixed_rows::<3>(0)) / (state.fixed_rows::<3>(0).norm().powi(3));
|
|
Vector6::new(state[3], state[4], state[5], acc[0], acc[1], acc[2])
|
|
}
|
|
let y0 = Vector6::new(
|
|
4.263868426884883e6,
|
|
5.146189057155391e6,
|
|
1.1310208421331816e6,
|
|
-5923.454461876975,
|
|
4496.802639690076,
|
|
1870.3893008991558,
|
|
);
|
|
|
|
// Integrate
|
|
let ode = ODE::new(&derivative, 0.0, 10.0 * period, y0, params);
|
|
let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-12);
|
|
let controller = PIController::new(0.37, 0.04, 10.0, 0.2, 1000.0, 0.9, 0.01);
|
|
let mut problem = Problem::new(ode, dp45, controller);
|
|
let solution = problem.solve();
|
|
|
|
assert_relative_eq!(
|
|
solution.times[solution.states.len() - 1],
|
|
10.0 * period,
|
|
max_relative = 1e-12
|
|
);
|
|
assert_relative_eq!(
|
|
solution.states[solution.states.len() - 1][0],
|
|
y0[0],
|
|
max_relative = 1e-9
|
|
);
|
|
assert_relative_eq!(
|
|
solution.states[solution.states.len() - 1][1],
|
|
y0[1],
|
|
max_relative = 1e-9
|
|
);
|
|
assert_relative_eq!(
|
|
solution.states[solution.states.len() - 1][2],
|
|
y0[2],
|
|
max_relative = 1e-9
|
|
);
|
|
assert_relative_eq!(
|
|
solution.states[solution.states.len() - 1][3],
|
|
y0[3],
|
|
max_relative = 1e-9
|
|
);
|
|
assert_relative_eq!(
|
|
solution.states[solution.states.len() - 1][4],
|
|
y0[4],
|
|
max_relative = 1e-9
|
|
);
|
|
assert_relative_eq!(
|
|
solution.states[solution.states.len() - 1][5],
|
|
y0[5],
|
|
max_relative = 1e-9
|
|
);
|
|
}
|
|
}
|