differential_equations/utils.rs
1use crate::{
2 Error,
3 traits::{CallBackData, Real, State},
4};
5
6/// Constrain the step size to be within the bounds of `h_min` and `h_max`.
7///
8/// If `h` is less than `h_min`, returns `h_min` with the same sign as `h`.
9/// If `h` is greater than `h_max`, returns `h_max` with the same sign as `h`.
10/// Otherwise, returns `h` unchanged.
11///
12/// # Arguments
13/// * `h` - Step size to constrain.
14/// * `h_min` - Minimum step size.
15/// * `h_max` - Maximum step size.
16///
17/// # Returns
18/// * The step size constrained to be within the bounds of `h_min` and `h_max`.
19///
20pub fn constrain_step_size<T: Real>(h: T, h_min: T, h_max: T) -> T {
21 // Determine the direction of the step size
22 let sign = h.signum();
23 // Bound the step size
24 if h.abs() < h_min {
25 sign * h_min
26 } else if h.abs() > h_max {
27 sign * h_max
28 } else {
29 h
30 }
31}
32
33/// Validate the step size parameters.
34///
35/// Checks the following:
36/// * `tf` cannot be equal to `t0`.
37/// * `h0` has the same sign as `tf - t0`.
38/// * `h_min` and `h_max` are non-negative.
39/// * `h_min` is less than or equal to `h_max`.
40/// * `|h0|` is greater than or equal to `h_min`.
41/// * `|h0|` is less than or equal to `h_max`.
42/// * `|h0|` is less than or equal to `|tf - t0|`.
43/// * `h0` is not zero.
44///
45/// If any of the checks fail, returns `Err(Error::BadInput)` with a descriptive message.
46/// Else returns `Ok(h0)` indicating the step size is valid.
47///
48/// # Arguments
49/// * `h0` - Initial step size.
50/// * `h_min` - Minimum step size.
51/// * `h_max` - Maximum step size.
52/// * `t0` - Initial time.
53/// * `tf` - Final time.
54///
55/// # Returns
56/// * `Result<Real, Error>` - Ok if all checks pass, Err if any check fails.
57///
58pub fn validate_step_size_parameters<T: Real, V: State<T>, D: CallBackData>(
59 h0: T,
60 h_min: T,
61 h_max: T,
62 t0: T,
63 tf: T,
64) -> Result<T, Error<T, V>> {
65 // Check if tf == t0
66 if tf == t0 {
67 return Err(Error::BadInput {
68 msg: format!(
69 "Invalid input: tf ({:?}) cannot be equal to t0 ({:?})",
70 tf, t0
71 ),
72 });
73 }
74
75 // Determine direction of the step size
76 let sign = (tf - t0).signum();
77
78 // Check h0 has same sign as tf - t0
79 if h0.signum() != sign {
80 return Err(Error::BadInput {
81 msg: format!(
82 "Invalid input: Initial step size ({:?}) must have the same sign as the integration direction (sign of tf - t0 = {:?})",
83 h0,
84 tf - t0
85 ),
86 });
87 }
88
89 // Check h_min and h_max bounds
90 if h_min < T::zero() {
91 return Err(Error::BadInput {
92 msg: format!(
93 "Invalid input: Minimum step size ({:?}) must be non-negative",
94 h_min
95 ),
96 });
97 }
98 if h_max < T::zero() {
99 return Err(Error::BadInput {
100 msg: format!(
101 "Invalid input: Maximum step size ({:?}) must be non-negative",
102 h_max
103 ),
104 });
105 }
106 if h_min > h_max {
107 return Err(Error::BadInput {
108 msg: format!(
109 "Invalid input: Minimum step size ({:?}) must be less than or equal to maximum step size ({:?})",
110 h_min, h_max
111 ),
112 });
113 }
114
115 // Check h0 bounds
116 if h0.abs() < h_min {
117 return Err(Error::BadInput {
118 msg: format!(
119 "Invalid input: Absolute value of initial step size ({:?}) must be greater than or equal to minimum step size ({:?})",
120 h0.abs(),
121 h_min
122 ),
123 });
124 }
125 if h0.abs() > h_max {
126 return Err(Error::BadInput {
127 msg: format!(
128 "Invalid input: Absolute value of initial step size ({:?}) must be less than or equal to maximum step size ({:?})",
129 h0.abs(),
130 h_max
131 ),
132 });
133 }
134
135 // Check h0 is not larger then integration interval
136 if h0.abs() > (tf - t0).abs() {
137 return Err(Error::BadInput {
138 msg: format!(
139 "Invalid input: Absolute value of initial step size ({:?}) must be less than or equal to the absolute value of the integration interval (tf - t0 = {:?})",
140 h0.abs(),
141 (tf - t0).abs()
142 ),
143 });
144 }
145
146 // Check h0 is not zero
147 if h0 == T::zero() {
148 return Err(Error::BadInput {
149 msg: format!("Invalid input: Initial step size ({:?}) cannot be zero", h0),
150 });
151 }
152
153 // Return Ok if all bounds are valid return the step size
154 Ok(h0)
155}