1use crate::{
2 error::Error,
3 interpolate::Interpolation,
4 linalg::{Matrix, lin_solve, lu_decomp},
5 methods::{Ordinary, h_init::InitialStepSize},
6 ode::{ODE, OrdinaryNumericalMethod},
7 stats::Evals,
8 status::Status,
9 traits::{Real, State},
10 utils::{constrain_step_size, validate_step_size_parameters},
11};
12
13use super::{BACKWARD_DIFFERENTIATION_FORMULA_ROWS, BackwardDifferentiationFormula, MAX_ORDER};
14
15impl<T: Real, Y: State<T>> BackwardDifferentiationFormula<Ordinary, T, Y> {
16 fn scalar(value: f64) -> T {
17 T::from_f64(value).expect("BackwardDifferentiationFormula constants must be representable as the solver scalar type")
18 }
19
20 fn order_scalar(order: usize) -> T {
21 T::from_usize(order)
22 .expect("BackwardDifferentiationFormula order constants must be representable as the solver scalar type")
23 }
24
25 fn init_coefficients(&mut self) {
26 let kappa = [
27 T::zero(),
28 Self::scalar(-0.1850),
29 -T::one() / Self::scalar(9.0),
30 Self::scalar(-0.0823),
31 Self::scalar(-0.0415),
32 T::zero(),
33 ];
34
35 self.gamma[0] = T::zero();
36 for order in 1..=MAX_ORDER {
37 self.gamma[order] = self.gamma[order - 1] + T::one() / Self::order_scalar(order);
38 self.alpha[order] = (T::one() - kappa[order]) * self.gamma[order];
39 self.error_const[order] =
40 kappa[order] * self.gamma[order] + T::one() / Self::order_scalar(order + 1);
41 }
42 }
43
44 fn weighted_rms_norm(&self, value: &Y, scale: &Y) -> T {
45 (value.component_div(scale).norm_squared() / Self::order_scalar(value.len())).sqrt()
46 }
47
48 fn scale_from(&self, y: &Y) -> Y {
49 let mut scale = y.zeros_like();
50 scale.map_components_mut(|i, s| {
51 *s = self.atol[i] + self.rtol[i] * y.get_component(i).abs();
52 });
53 scale
54 }
55
56 fn predict(&self, order: usize) -> Y {
57 let mut y_predict = self.y.zeros_like();
58 for i in 0..=order {
59 y_predict.add_scaled(T::one(), &self.d[i]);
60 }
61 y_predict
62 }
63
64 fn compute_psi(&self, order: usize) -> Y {
65 let mut psi = self.y.zeros_like();
66 for i in 1..=order {
67 psi.add_scaled(self.gamma[i], &self.d[i]);
68 }
69 psi.scaled(T::one() / self.alpha[order])
70 }
71
72 fn compute_r(
73 order: usize,
74 factor: T,
75 ) -> [[T; BACKWARD_DIFFERENTIATION_FORMULA_ROWS]; BACKWARD_DIFFERENTIATION_FORMULA_ROWS] {
76 let mut r = [[T::zero(); BACKWARD_DIFFERENTIATION_FORMULA_ROWS];
77 BACKWARD_DIFFERENTIATION_FORMULA_ROWS];
78 for j in 0..=order {
79 r[0][j] = T::one();
80 }
81
82 for i in 1..=order {
83 for j in 1..=order {
84 r[i][j] = (Self::order_scalar(i - 1) - factor * Self::order_scalar(j))
85 / Self::order_scalar(i);
86 }
87 }
88
89 for i in 1..=order {
90 for j in 1..=order {
91 r[i][j] *= r[i - 1][j];
92 }
93 }
94
95 r
96 }
97
98 fn change_d(&mut self, order: usize, factor: T) {
99 if factor == T::one() {
100 return;
101 }
102
103 let r = Self::compute_r(order, factor);
104 let u = Self::compute_r(order, T::one());
105 let mut ru = [[T::zero(); BACKWARD_DIFFERENTIATION_FORMULA_ROWS];
106 BACKWARD_DIFFERENTIATION_FORMULA_ROWS];
107
108 for i in 0..=order {
109 for j in 0..=order {
110 let mut sum = T::zero();
111 for (k, u_row) in u.iter().enumerate().take(order + 1) {
112 sum += r[i][k] * u_row[j];
113 }
114 ru[i][j] = sum;
115 }
116 }
117
118 let old = self.d.clone();
119 for i in 0..=order {
120 let mut transformed = self.y.zeros_like();
121 for j in 0..=order {
122 transformed.add_scaled(ru[j][i], &old[j]);
123 }
124 self.d[i] = transformed;
125 }
126 }
127
128 fn factor_matrix(&mut self, c: T) -> bool {
129 let dim = self.y.len();
130 for i in 0..dim {
131 for j in 0..dim {
132 self.newton_matrix[(i, j)] = if i == j {
133 T::one() - c * self.jacobian[(i, j)]
134 } else {
135 -c * self.jacobian[(i, j)]
136 };
137 }
138 }
139
140 lu_decomp(&mut self.newton_matrix, &mut self.ip).is_ok()
141 }
142
143 fn solve_system<F>(
144 &mut self,
145 ode: &F,
146 t_new: T,
147 y_predict: Y,
148 c: T,
149 psi: Y,
150 scale: Y,
151 evals: &mut Evals,
152 ) -> (bool, usize, Y, Y)
153 where
154 F: ODE<T, Y> + ?Sized,
155 {
156 let mut d = self.y.zeros_like();
157 let mut y = y_predict;
158 let mut dy_norm_old: Option<T> = None;
159
160 for iter in 0..self.max_newton_iter {
161 ode.diff(t_new, &y, &mut self.dydt);
162 evals.function += 1;
163
164 let mut dy = self.dydt.scaled(c);
165 dy.add_scaled(-T::one(), &psi);
166 dy.add_scaled(-T::one(), &d);
167 lin_solve(&self.newton_matrix, &mut dy, &self.ip);
168 let dy_norm = self.weighted_rms_norm(&dy, &scale);
169
170 let rate = dy_norm_old.map(|old| dy_norm / old);
171 if let Some(rate) = rate {
172 let remaining = self.max_newton_iter - iter;
173 if rate >= T::one()
174 || rate.powf(Self::order_scalar(remaining)) / (T::one() - rate) * dy_norm
175 > self.newton_tol
176 {
177 return (false, iter + 1, y, d);
178 }
179 }
180
181 y.add_scaled(T::one(), &dy);
182 d.add_scaled(T::one(), &dy);
183
184 if dy_norm == T::zero()
185 || rate.is_some_and(|rate| rate / (T::one() - rate) * dy_norm < self.newton_tol)
186 {
187 return (true, iter + 1, y, d);
188 }
189
190 dy_norm_old = Some(dy_norm);
191 }
192
193 (false, self.max_newton_iter, y, d)
194 }
195
196 fn step_too_small(&self) -> bool {
197 let min_step = self
198 .h_min
199 .max(T::default_epsilon() * Self::scalar(10.0) * (self.t.abs() + T::one()));
200 self.h.abs() < min_step
201 }
202
203 fn apply_step_factor(&mut self, order: usize, factor: T) {
204 self.h = constrain_step_size(self.h * factor, self.h_min, self.h_max);
205 self.h = (self.filter)(self.h);
206 self.change_d(order, factor);
207 self.n_equal_steps = 0;
208 self.lu_valid = false;
209 }
210}
211
212impl<T: Real, Y: State<T>> OrdinaryNumericalMethod<T, Y>
213 for BackwardDifferentiationFormula<Ordinary, T, Y>
214{
215 fn init<F>(&mut self, ode: &F, t0: T, tf: T, y0: &Y) -> Result<Evals, Error<T, Y>>
216 where
217 F: ODE<T, Y> + ?Sized,
218 {
219 let mut evals = Evals::new();
220
221 if self.h0 == T::zero() {
222 self.h0 = InitialStepSize::<Ordinary>::compute(
223 ode, t0, tf, y0, 1, &self.rtol, &self.atol, self.h_min, self.h_max, &mut evals,
224 );
225 }
226
227 self.h = validate_step_size_parameters(self.h0, self.h_min, self.h_max, t0, tf)?;
228 self.h = (self.filter)(self.h);
229 self.init_coefficients();
230 if self.newton_tol == T::zero() {
231 let rtol = self.rtol.average();
232 self.newton_tol = (Self::scalar(10.0) * T::default_epsilon() / rtol)
233 .max(Self::scalar(0.03).min(rtol.sqrt()));
234 }
235
236 self.t = t0;
237 self.y = y0.clone();
238 self.t_prev = t0;
239 self.y_prev = y0.clone();
240 self.h_prev = self.h;
241 self.order = 1;
242 self.n_equal_steps = 0;
243 self.steps = 0;
244 self.lu_valid = false;
245
246 ode.diff(t0, y0, &mut self.dydt);
247 evals.function += 1;
248
249 let dim = y0.len();
250 self.d = core::array::from_fn(|_| y0.zeros_like());
251 self.d[0] = y0.clone();
252 self.d[1] = self.dydt.scaled(self.h);
253 self.jacobian = Matrix::zeros(dim, dim);
254 self.newton_matrix = Matrix::zeros(dim, dim);
255 self.ip = vec![0; dim];
256
257 ode.jacobian(t0, y0, &mut self.jacobian);
258 evals.jacobian += 1;
259
260 self.status = Status::Initialized;
261 Ok(evals)
262 }
263
264 fn step<F>(&mut self, ode: &F) -> Result<Evals, Error<T, Y>>
265 where
266 F: ODE<T, Y> + ?Sized,
267 {
268 let mut evals = Evals::new();
269
270 if self.step_too_small() {
271 let e = Error::StepSize {
272 t: self.t,
273 y: self.y.clone(),
274 };
275 self.status = Status::Error(e.clone());
276 return Err(e);
277 }
278
279 if self.steps >= self.max_steps {
280 let e = Error::MaxSteps {
281 t: self.t,
282 y: self.y.clone(),
283 };
284 self.status = Status::Error(e.clone());
285 return Err(e);
286 }
287
288 let mut order = self.order;
289 let mut step_accepted = false;
290 let mut accepted_y = self.y.clone();
291 let mut accepted_d = self.y.zeros_like();
292 let mut accepted_error_norm = T::zero();
293 let mut accepted_safety = T::one();
294
295 while !step_accepted {
296 if self.step_too_small() {
297 let e = Error::StepSize {
298 t: self.t,
299 y: self.y.clone(),
300 };
301 self.status = Status::Error(e.clone());
302 return Err(e);
303 }
304
305 let t_new = self.t + self.h;
306 let y_predict = self.predict(order);
307 let scale_predict = self.scale_from(&y_predict);
308 let psi = self.compute_psi(order);
309 let c = self.h / self.alpha[order];
310
311 let mut converged = false;
312 let mut n_iter = self.max_newton_iter;
313 let mut y_new = y_predict.clone();
314 let mut d = self.y.zeros_like();
315 let mut jacobian_current = false;
316
317 while !converged {
318 if !self.lu_valid {
319 if !self.factor_matrix(c) {
320 converged = false;
321 break;
322 }
323 self.lu_valid = true;
324 }
325
326 (converged, n_iter, y_new, d) = self.solve_system(
327 ode,
328 t_new,
329 y_predict.clone(),
330 c,
331 psi.clone(),
332 scale_predict.clone(),
333 &mut evals,
334 );
335
336 if !converged {
337 if jacobian_current {
338 break;
339 }
340 ode.jacobian(t_new, &y_predict, &mut self.jacobian);
341 evals.jacobian += 1;
342 self.lu_valid = false;
343 jacobian_current = true;
344 }
345 }
346
347 if !converged {
348 self.apply_step_factor(order, Self::scalar(0.5));
349 continue;
350 }
351
352 let safety = Self::scalar(0.9) * Self::order_scalar(2 * self.max_newton_iter + 1)
353 / Self::order_scalar(2 * self.max_newton_iter + n_iter);
354 let scale = self.scale_from(&y_new);
355 let error = d.scaled(self.error_const[order]);
356 let error_norm = self.weighted_rms_norm(&error, &scale);
357
358 if error_norm > T::one() {
359 let factor = self
360 .min_scale
361 .max(safety * error_norm.powf(-T::one() / Self::order_scalar(order + 1)));
362 self.apply_step_factor(order, factor);
363 } else {
364 step_accepted = true;
365 accepted_y = y_new;
366 accepted_d = d;
367 accepted_error_norm = error_norm;
368 accepted_safety = safety;
369 }
370 }
371
372 self.t_prev = self.t;
373 self.y_prev = self.y.clone();
374 self.h_prev = self.h;
375 self.t += self.h;
376 self.y = accepted_y;
377 self.steps += 1;
378 self.n_equal_steps += 1;
379
380 self.d[order + 2] = accepted_d.minus(&self.d[order + 1]);
381 self.d[order + 1] = accepted_d.clone();
382 for i in (0..=order).rev() {
383 let next = self.d[i + 1].clone();
384 self.d[i].add_scaled(T::one(), &next);
385 }
386
387 if self.n_equal_steps > order {
388 let scale = self.scale_from(&self.y);
389 let error_m_norm = if order > 1 {
390 self.weighted_rms_norm(&self.d[order].scaled(self.error_const[order - 1]), &scale)
391 } else {
392 T::infinity()
393 };
394 let error_p_norm = if order < self.max_order {
395 self.weighted_rms_norm(
396 &self.d[order + 2].scaled(self.error_const[order + 1]),
397 &scale,
398 )
399 } else {
400 T::infinity()
401 };
402
403 let current_factor =
404 accepted_error_norm.powf(-T::one() / Self::order_scalar(order + 1));
405 let lower_factor = error_m_norm.powf(-T::one() / Self::order_scalar(order));
406 let higher_factor = error_p_norm.powf(-T::one() / Self::order_scalar(order + 2));
407
408 let mut best_factor = current_factor;
409 let mut delta_order: isize = 0;
410 if lower_factor > best_factor {
411 best_factor = lower_factor;
412 delta_order = -1;
413 }
414 if higher_factor > best_factor {
415 best_factor = higher_factor;
416 delta_order = 1;
417 }
418
419 order = (order as isize + delta_order) as usize;
420 self.order = order;
421
422 let factor = self.max_scale.min(accepted_safety * best_factor);
423 self.apply_step_factor(order, factor);
424 }
425
426 self.status = Status::Solving;
427 Ok(evals)
428 }
429
430 fn t(&self) -> T {
431 self.t
432 }
433
434 fn y(&self) -> &Y {
435 &self.y
436 }
437
438 fn t_prev(&self) -> T {
439 self.t_prev
440 }
441
442 fn y_prev(&self) -> &Y {
443 &self.y_prev
444 }
445
446 fn h(&self) -> T {
447 self.h
448 }
449
450 fn set_h(&mut self, h: T) {
451 let h = (self.filter)(h);
452 if self.h != T::zero() {
453 let factor = h / self.h;
454 self.change_d(self.order, factor);
455 self.n_equal_steps = 0;
456 self.lu_valid = false;
457 }
458 self.h = h;
459 }
460
461 fn status(&self) -> &Status<T, Y> {
462 &self.status
463 }
464
465 fn set_status(&mut self, status: Status<T, Y>) {
466 self.status = status;
467 }
468}
469
470impl<T: Real, Y: State<T>> Interpolation<T, Y> for BackwardDifferentiationFormula<Ordinary, T, Y> {
471 fn interpolate(&mut self, t_interp: T) -> Result<Y, Error<T, Y>> {
472 if t_interp < self.t_prev || t_interp > self.t {
473 return Err(Error::OutOfBounds {
474 t_interp,
475 t_prev: self.t_prev,
476 t_curr: self.t,
477 });
478 }
479
480 if self.h_prev == T::zero() {
481 return Ok(self.y.clone());
482 }
483
484 let mut y_interp = self.d[0].clone();
485 let mut product = T::one();
486 for j in 1..=self.order {
487 let shift = self.t - self.h_prev * Self::order_scalar(j - 1);
488 let denom = self.h_prev * Self::order_scalar(j);
489 product = product * (t_interp - shift) / denom;
490 y_interp.add_scaled(product, &self.d[j]);
491 }
492
493 Ok(y_interp)
494 }
495}
496
497#[cfg(all(test, feature = "nalgebra"))]
498mod tests {
499 use super::*;
500 use crate::prelude::*;
501 use nalgebra::{SVector, vector};
502
503 struct ExponentialDecay {
504 k: f64,
505 }
506
507 impl ODE<f64, SVector<f64, 1>> for ExponentialDecay {
508 fn diff(&self, _t: f64, y: &SVector<f64, 1>, dydt: &mut SVector<f64, 1>) {
509 dydt[0] = -self.k * y[0];
510 }
511
512 fn jacobian(&self, _t: f64, _y: &SVector<f64, 1>, j: &mut Matrix<f64>) {
513 j[(0, 0)] = -self.k;
514 }
515 }
516
517 #[test]
518 fn backward_differentiation_formula_tracks_exponential_decay() {
519 let system = ExponentialDecay { k: 1.0 };
520 let backward_differentiation_formula: BackwardDifferentiationFormula<
521 Ordinary,
522 f64,
523 SVector<f64, 1>,
524 > = BackwardDifferentiationFormula::adaptive()
525 .rtol(1e-7)
526 .atol(1e-9);
527
528 let results = IVP::ode(&system, 0.0, 1.0, vector![1.0])
529 .method(backward_differentiation_formula)
530 .solve()
531 .unwrap();
532
533 let actual = results.y.last().unwrap()[0];
534 let expected = (-1.0_f64).exp();
535 assert!((actual - expected).abs() < 1e-6);
536 }
537
538 #[test]
539 fn backward_differentiation_formula_can_run_at_order_one() {
540 let system = ExponentialDecay { k: 1.0 };
541 let backward_differentiation_formula: BackwardDifferentiationFormula<
542 Ordinary,
543 f64,
544 SVector<f64, 1>,
545 > = BackwardDifferentiationFormula::adaptive()
546 .max_order(1)
547 .rtol(1e-7)
548 .atol(1e-9);
549
550 let results = IVP::ode(&system, 0.0, 1.0, vector![1.0])
551 .method(backward_differentiation_formula)
552 .solve()
553 .unwrap();
554
555 let actual = results.y.last().unwrap()[0];
556 let expected = (-1.0_f64).exp();
557 assert!(
558 (actual - expected).abs() < 1e-4,
559 "got {actual}, expected {expected}"
560 );
561 }
562}