1#[allow(clippy::excessive_precision, clippy::unreadable_literal)]
5mod coefficients {
6 pub const C: [f64; 16] = [
7 0.0,
8 0.526001519587677318785587544488e-01,
9 0.789002279381515978178381316732e-01,
10 0.118350341907227396726757197510e+00,
11 0.281649658092772603273242802490e+00,
12 0.333333333333333333333333333333e+00,
13 0.25e+00,
14 0.307692307692307692307692307692e+00,
15 0.651282051282051282051282051282e+00,
16 0.6e+00,
17 0.857142857142857142857142857142e+00,
18 1.0,
19 0.0,
20 0.1e+00,
21 0.2e+00,
22 0.777777777777777777777777777778e+00,
23 ];
24 pub const A21: f64 = 5.26001519587677318785587544488e-2;
25 pub const A31: f64 = 1.97250569845378994544595329183e-2;
26 pub const A32: f64 = 5.91751709536136983633785987549e-2;
27 pub const A41: f64 = 2.95875854768068491816892993775e-2;
28 pub const A43: f64 = 8.87627564304205475450678981324e-2;
29 pub const A51: f64 = 2.41365134159266685502369798665e-1;
30 pub const A53: f64 = -8.84549479328286085344864962717e-1;
31 pub const A54: f64 = 9.24834003261792003115737966543e-1;
32 pub const A61: f64 = 3.7037037037037037037037037037e-2;
33 pub const A64: f64 = 1.70828608729473871279604482173e-1;
34 pub const A65: f64 = 1.25467687566822425016691814123e-1;
35 pub const A71: f64 = 3.7109375e-2;
36 pub const A74: f64 = 1.70252211019544039314978060272e-1;
37 pub const A75: f64 = 6.02165389804559606850219397283e-2;
38 pub const A76: f64 = -1.7578125e-2;
39 pub const A81: f64 = 3.70920001185047927108779319836e-2;
40 pub const A84: f64 = 1.70383925712239993810214054705e-1;
41 pub const A85: f64 = 1.07262030446373284651809199168e-1;
42 pub const A86: f64 = -1.53194377486244017527936158236e-2;
43 pub const A87: f64 = 8.27378916381402288758473766002e-3;
44 pub const A91: f64 = 6.24110958716075717114429577812e-1;
45 pub const A94: f64 = -3.36089262944694129406857109825e+0;
46 pub const A95: f64 = -8.68219346841726006818189891453e-1;
47 pub const A96: f64 = 2.75920996994467083049415600797e+1;
48 pub const A97: f64 = 2.01540675504778934086186788979e+1;
49 pub const A98: f64 = -4.34898841810699588477366255144e+1;
50 pub const A101: f64 = 4.77662536438264365890433908527e-1;
51 pub const A104: f64 = -2.48811461997166764192642586468e+0;
52 pub const A105: f64 = -5.90290826836842996371446475743e-1;
53 pub const A106: f64 = 2.12300514481811942347288949897e+1;
54 pub const A107: f64 = 1.52792336328824235832596922938e+1;
55 pub const A108: f64 = -3.32882109689848629194453265587e+1;
56 pub const A109: f64 = -2.03312017085086261358222928593e-2;
57 pub const A111: f64 = -9.3714243008598732571704021658e-1;
58 pub const A114: f64 = 5.18637242884406370830023853209e+0;
59 pub const A115: f64 = 1.09143734899672957818500254654e+0;
60 pub const A116: f64 = -8.14978701074692612513997267357e+0;
61 pub const A117: f64 = -1.85200656599969598641566180701e+1;
62 pub const A118: f64 = 2.27394870993505042818970056734e+1;
63 pub const A119: f64 = 2.49360555267965238987089396762e+0;
64 pub const A1110: f64 = -3.0467644718982195003823669022e+0;
65 pub const A121: f64 = 2.27331014751653820792359768449e+0;
66 pub const A124: f64 = -1.05344954667372501984066689879e+1;
67 pub const A125: f64 = -2.00087205822486249909675718444e+0;
68 pub const A126: f64 = -1.79589318631187989172765950534e+1;
69 pub const A127: f64 = 2.79488845294199600508499808837e+1;
70 pub const A128: f64 = -2.85899827713502369474065508674e+0;
71 pub const A129: f64 = -8.87285693353062954433549289258e+0;
72 pub const A1210: f64 = 1.23605671757943030647266201528e+1;
73 pub const A1211: f64 = 6.43392746015763530355970484046e-1;
74 pub const B: [f64; 12] = [
75 5.42937341165687622380535766363e-2,
76 0.0,
77 0.0,
78 0.0,
79 0.0,
80 4.45031289275240888144113950566e+0,
81 1.89151789931450038304281599044e+0,
82 -5.8012039600105847814672114227e+0,
83 3.1116436695781989440891606237e-1,
84 -1.52160949662516078556178806805e-1,
85 2.01365400804030348374776537501e-1,
86 4.47106157277725905176885569043e-2,
87 ];
88 pub const ER: [f64; 12] = [
89 0.1312004499419488073250102996e-1,
90 0.0,
91 0.0,
92 0.0,
93 0.0,
94 -0.1225156446376204440720569753e+1,
95 -0.4957589496572501915214079952e+0,
96 0.1664377182454986536961530415e+1,
97 -0.3503288487499736816886487290e+0,
98 0.3341791187130174790297318841e+0,
99 0.8192320648511571246570742613e-1,
100 -0.2235530786388629525884427845e-1,
101 ];
102 pub const BHH: [f64; 3] = [
103 0.244094488188976377952755905512e+0,
104 0.733846688281611857341361741547e+0,
105 0.220588235294117647058823529412e-1,
106 ];
107}
108
109use coefficients::*;
110
111#[derive(Debug, Clone)]
112pub struct Dop853Options {
113 pub rtol: f64,
114 pub atol: f64,
115 pub max_step: f64,
116 pub initial_step: f64,
117}
118
119impl Default for Dop853Options {
120 fn default() -> Self {
121 Self {
122 rtol: 1e-10,
123 atol: 1e-13,
124 max_step: 5e-4,
125 initial_step: 0.0,
126 }
127 }
128}
129
130#[derive(Debug, Clone)]
131pub struct Dop853Result {
132 pub t: Vec<f64>,
133 pub y: Vec<Vec<f64>>,
134 pub success: bool,
135 pub n_steps: usize,
136 pub n_evals: usize,
137}
138
139pub fn integrate_dop853<F>(
140 rhs: F,
141 t_span: (f64, f64),
142 y0: &[f64],
143 t_eval: &[f64],
144 opts: &Dop853Options,
145) -> Dop853Result
146where
147 F: Fn(f64, &[f64], &mut [f64]),
148{
149 let n = y0.len();
150 let t0 = t_span.0;
151 let t_end = t_span.1;
152 let mut t = t0;
153 let mut y = y0.to_vec();
154 let mut h = if opts.initial_step > 0.0 {
155 opts.initial_step.min(opts.max_step)
156 } else {
157 initial_step_size(&rhs, t0, &y, opts)
158 };
159
160 let mut result_t: Vec<f64> = Vec::with_capacity(t_eval.len());
161 let mut result_y: Vec<Vec<f64>> = Vec::with_capacity(t_eval.len());
162 let mut eval_idx = 0;
163 let mut n_steps = 0usize;
164 let mut n_evals = 0usize;
165 let mut step_rejected = false;
166
167 while eval_idx < t_eval.len() && t_eval[eval_idx] <= t0 + 1e-15 * t0.abs().max(1.0) {
168 result_t.push(t_eval[eval_idx]);
169 result_y.push(y.clone());
170 eval_idx += 1;
171 }
172
173 let mut k = vec![vec![0.0; n]; 13];
174 rhs(t, &y, &mut k[0]);
175 n_evals += 1;
176 let max_steps = 10_000_000;
177
178 while t < t_end && n_steps < max_steps {
179 if t + h > t_end {
180 h = t_end - t;
181 }
182 if h < 1e-15 * t.abs().max(1.0) {
183 break;
184 }
185
186 let mut ys = vec![0.0; n];
187 for i in 0..n {
188 ys[i] = y[i] + h * A21 * k[0][i];
189 }
190 rhs(t + C[1] * h, &ys, &mut k[1]);
191 for i in 0..n {
192 ys[i] = y[i] + h * (A31 * k[0][i] + A32 * k[1][i]);
193 }
194 rhs(t + C[2] * h, &ys, &mut k[2]);
195 for i in 0..n {
196 ys[i] = y[i] + h * (A41 * k[0][i] + A43 * k[2][i]);
197 }
198 rhs(t + C[3] * h, &ys, &mut k[3]);
199 for i in 0..n {
200 ys[i] = y[i] + h * (A51 * k[0][i] + A53 * k[2][i] + A54 * k[3][i]);
201 }
202 rhs(t + C[4] * h, &ys, &mut k[4]);
203 for i in 0..n {
204 ys[i] = y[i] + h * (A61 * k[0][i] + A64 * k[3][i] + A65 * k[4][i]);
205 }
206 rhs(t + C[5] * h, &ys, &mut k[5]);
207 for i in 0..n {
208 ys[i] = y[i] + h * (A71 * k[0][i] + A74 * k[3][i] + A75 * k[4][i] + A76 * k[5][i]);
209 }
210 rhs(t + C[6] * h, &ys, &mut k[6]);
211 for i in 0..n {
212 ys[i] = y[i]
213 + h * (A81 * k[0][i]
214 + A84 * k[3][i]
215 + A85 * k[4][i]
216 + A86 * k[5][i]
217 + A87 * k[6][i]);
218 }
219 rhs(t + C[7] * h, &ys, &mut k[7]);
220 for i in 0..n {
221 ys[i] = y[i]
222 + h * (A91 * k[0][i]
223 + A94 * k[3][i]
224 + A95 * k[4][i]
225 + A96 * k[5][i]
226 + A97 * k[6][i]
227 + A98 * k[7][i]);
228 }
229 rhs(t + C[8] * h, &ys, &mut k[8]);
230 for i in 0..n {
231 ys[i] = y[i]
232 + h * (A101 * k[0][i]
233 + A104 * k[3][i]
234 + A105 * k[4][i]
235 + A106 * k[5][i]
236 + A107 * k[6][i]
237 + A108 * k[7][i]
238 + A109 * k[8][i]);
239 }
240 rhs(t + C[9] * h, &ys, &mut k[9]);
241 for i in 0..n {
242 ys[i] = y[i]
243 + h * (A111 * k[0][i]
244 + A114 * k[3][i]
245 + A115 * k[4][i]
246 + A116 * k[5][i]
247 + A117 * k[6][i]
248 + A118 * k[7][i]
249 + A119 * k[8][i]
250 + A1110 * k[9][i]);
251 }
252 rhs(t + C[10] * h, &ys, &mut k[10]);
253 for i in 0..n {
254 ys[i] = y[i]
255 + h * (A121 * k[0][i]
256 + A124 * k[3][i]
257 + A125 * k[4][i]
258 + A126 * k[5][i]
259 + A127 * k[6][i]
260 + A128 * k[7][i]
261 + A129 * k[8][i]
262 + A1210 * k[9][i]
263 + A1211 * k[10][i]);
264 }
265 rhs(t + h, &ys, &mut k[11]);
266 n_evals += 11;
267
268 let mut y_new = vec![0.0; n];
269 for i in 0..n {
270 y_new[i] = y[i]
271 + h * (B[0] * k[0][i]
272 + B[5] * k[5][i]
273 + B[6] * k[6][i]
274 + B[7] * k[7][i]
275 + B[8] * k[8][i]
276 + B[9] * k[9][i]
277 + B[10] * k[10][i]
278 + B[11] * k[11][i]);
279 }
280
281 rhs(t + h, &y_new, &mut k[12]);
282 n_evals += 1;
283
284 let mut err_sq = 0.0;
285 let mut err2_sq = 0.0;
286 for i in 0..n {
287 let sk = opts.atol + opts.rtol * y[i].abs().max(y_new[i].abs());
288 let er_i = ER[0] * k[0][i]
289 + ER[5] * k[5][i]
290 + ER[6] * k[6][i]
291 + ER[7] * k[7][i]
292 + ER[8] * k[8][i]
293 + ER[9] * k[9][i]
294 + ER[10] * k[10][i]
295 + ER[11] * k[11][i];
296 let slope_i = (y_new[i] - y[i]) / h;
297 let bhh_i = slope_i - BHH[0] * k[0][i] - BHH[1] * k[8][i] - BHH[2] * k[11][i];
298 err_sq += (er_i / sk) * (er_i / sk);
299 err2_sq += (bhh_i / sk) * (bhh_i / sk);
300 }
301
302 let denom = err_sq + 0.01 * err2_sq;
303 let err = if denom > 0.0 {
304 h.abs() * err_sq / (denom * n as f64).sqrt()
305 } else {
306 0.0
307 };
308
309 if err <= 1.0 {
310 let t_new = t + h;
311 n_steps += 1;
312 while eval_idx < t_eval.len()
313 && t_eval[eval_idx] <= t_new + 1e-15 * t_new.abs().max(1.0)
314 {
315 let te = t_eval[eval_idx];
316 if (te - t_new).abs() < 1e-15 * t_new.abs().max(1.0) {
317 result_t.push(te);
318 result_y.push(y_new.clone());
319 } else {
320 let theta = (te - t) / h;
321 let mut y_interp = vec![0.0; n];
322 for i in 0..n {
323 y_interp[i] = y[i] * (1.0 - theta)
324 + y_new[i] * theta
325 + theta
326 * (theta - 1.0)
327 * ((1.0 - 2.0 * theta) * (y_new[i] - y[i])
328 + (theta - 1.0) * h * k[0][i]
329 + theta * h * k[12][i]);
330 }
331 result_t.push(te);
332 result_y.push(y_interp);
333 }
334 eval_idx += 1;
335 }
336
337 t = t_new;
338 y = y_new;
339 k[0] = k[12].clone();
340 let fac = if err > 0.0 {
341 0.9 * err.powf(-1.0 / 8.0)
342 } else {
343 5.0
344 };
345 let fac = if step_rejected {
346 fac.min(1.0)
347 } else {
348 fac.clamp(0.2, 10.0)
349 };
350 h *= fac;
351 h = h.min(opts.max_step);
352 step_rejected = false;
353 } else {
354 let fac = 0.9 * err.powf(-1.0 / 8.0);
355 h *= fac.max(0.2);
356 step_rejected = true;
357 }
358 }
359
360 Dop853Result {
361 t: result_t,
362 y: result_y,
363 success: n_steps < max_steps && eval_idx == t_eval.len(),
364 n_steps,
365 n_evals,
366 }
367}
368
369fn initial_step_size<F>(rhs: &F, t0: f64, y0: &[f64], opts: &Dop853Options) -> f64
370where
371 F: Fn(f64, &[f64], &mut [f64]),
372{
373 let n = y0.len();
374 let mut f0 = vec![0.0; n];
375 rhs(t0, y0, &mut f0);
376 let mut d0 = 0.0;
377 let mut d1 = 0.0;
378 for i in 0..n {
379 let sk = opts.atol + opts.rtol * y0[i].abs();
380 d0 += (y0[i] / sk) * (y0[i] / sk);
381 d1 += (f0[i] / sk) * (f0[i] / sk);
382 }
383 d0 = (d0 / n as f64).sqrt();
384 d1 = (d1 / n as f64).sqrt();
385 let h0 = if d0 < 1e-5 || d1 < 1e-5 {
386 1e-6
387 } else {
388 0.01 * d0 / d1
389 };
390 let h0 = h0.min(opts.max_step);
391
392 let mut y1 = vec![0.0; n];
393 for i in 0..n {
394 y1[i] = y0[i] + h0 * f0[i];
395 }
396 let mut f1 = vec![0.0; n];
397 rhs(t0 + h0, &y1, &mut f1);
398 let mut d2 = 0.0;
399 for i in 0..n {
400 let sk = opts.atol + opts.rtol * y0[i].abs();
401 d2 += ((f1[i] - f0[i]) / sk) * ((f1[i] - f0[i]) / sk);
402 }
403 d2 = (d2 / n as f64).sqrt() / h0;
404
405 let h1 = if d1.max(d2) <= 1e-15 {
406 (h0 * 1e-3).max(1e-6)
407 } else {
408 (0.01 / d1.max(d2)).powf(1.0 / 8.0)
409 };
410 h0.min(100.0 * h1).min(opts.max_step)
411}
412
413#[cfg(test)]
414mod tests {
415 use super::*;
416
417 #[test]
418 fn test_harmonic_oscillator() {
419 let rhs = |_t: f64, y: &[f64], dydt: &mut [f64]| {
420 dydt[0] = y[1];
421 dydt[1] = -y[0];
422 };
423 let y0 = [1.0, 0.0];
424 let t_end = 10.0;
425 let n_eval = 1000;
426 let t_eval: Vec<f64> = (0..=n_eval)
427 .map(|i| i as f64 * t_end / n_eval as f64)
428 .collect();
429 let result = integrate_dop853(
430 rhs,
431 (0.0, t_end),
432 &y0,
433 &t_eval,
434 &Dop853Options {
435 rtol: 1e-10,
436 atol: 1e-13,
437 max_step: 0.05,
438 initial_step: 0.0,
439 },
440 );
441 assert!(result.success, "DOP853 failed: steps={}", result.n_steps);
442 assert_eq!(result.t.len(), t_eval.len());
443 for (i, &ti) in result.t.iter().enumerate() {
444 let exact_y = ti.cos();
445 let err_y = (result.y[i][0] - exact_y).abs();
446 assert!(err_y < 1e-7, "t={ti:.4}: y err = {err_y:.2e}");
447 }
448 let final_y = result.y.last().unwrap();
449 let exact_final = t_end.cos();
450 let err = (final_y[0] - exact_final).abs();
451 assert!(err < 1e-9, "final err = {err:.2e}");
452 }
453
454 #[test]
455 fn test_exponential_decay() {
456 let rhs = |_t: f64, y: &[f64], dydt: &mut [f64]| {
457 dydt[0] = -y[0];
458 };
459 let y0 = [1.0];
460 let t_end = 5.0;
461 let t_eval: Vec<f64> = (0..=50).map(|i| i as f64 * t_end / 50.0).collect();
462 let result = integrate_dop853(rhs, (0.0, t_end), &y0, &t_eval, &Dop853Options::default());
463 assert!(result.success, "DOP853 failed: steps={}", result.n_steps);
464 for (i, &ti) in result.t.iter().enumerate() {
465 let exact = (-ti).exp();
466 let err = (result.y[i][0] - exact).abs();
467 assert!(err < 1e-10, "t={ti:.2}: err = {err:.2e}");
468 }
469 }
470
471 #[test]
472 fn test_energy_conservation() {
473 let rhs = |_t: f64, y: &[f64], dydt: &mut [f64]| {
474 dydt[0] = -y[1];
475 dydt[1] = y[0];
476 };
477 let y0 = [1.0, 0.0];
478 let t_end = 100.0;
479 let t_eval: Vec<f64> = (0..=1000).map(|i| i as f64 * t_end / 1000.0).collect();
480 let result = integrate_dop853(
481 rhs,
482 (0.0, t_end),
483 &y0,
484 &t_eval,
485 &Dop853Options {
486 rtol: 1e-12,
487 atol: 1e-14,
488 max_step: 0.01,
489 initial_step: 0.0,
490 },
491 );
492 assert!(result.success, "DOP853 failed: steps={}", result.n_steps);
493 let e0 = y0[0] * y0[0] + y0[1] * y0[1];
494 for y in &result.y {
495 let e = y[0] * y[0] + y[1] * y[1];
496 let err = (e - e0).abs();
497 assert!(err < 1e-10, "energy err = {err:.2e}");
498 }
499 }
500
501 #[test]
502 fn test_logistic_growth_regression() {
503 let rhs = |_t: f64, y: &[f64], dydt: &mut [f64]| {
504 let r = 1.5;
505 let k = 2.0;
506 dydt[0] = r * y[0] * (1.0 - y[0] / k);
507 };
508 let y0 = [0.25];
509 let t_end = 3.0;
510 let t_eval: Vec<f64> = (0..=30).map(|i| i as f64 * t_end / 30.0).collect();
511 let result = integrate_dop853(
512 rhs,
513 (0.0, t_end),
514 &y0,
515 &t_eval,
516 &Dop853Options {
517 rtol: 1e-11,
518 atol: 1e-13,
519 max_step: 0.05,
520 initial_step: 0.0,
521 },
522 );
523 assert!(result.success, "DOP853 failed: steps={}", result.n_steps);
524 assert_eq!(result.t, t_eval);
525 for (i, &ti) in result.t.iter().enumerate() {
526 let exact = 2.0 / (1.0 + 7.0 * (-1.5 * ti).exp());
527 let err = (result.y[i][0] - exact).abs();
528 assert!(err < 1e-8, "t={ti:.3}: err = {err:.2e}");
529 }
530 }
531
532 #[test]
533 fn test_includes_initial_state_when_t_eval_starts_at_t0() {
534 let rhs = |_t: f64, y: &[f64], dydt: &mut [f64]| {
535 dydt[0] = -2.0 * y[0];
536 };
537 let t_eval = [0.0, 0.25, 0.5];
538 let y0 = [3.0];
539 let result = integrate_dop853(rhs, (0.0, 0.5), &y0, &t_eval, &Dop853Options::default());
540 assert!(result.success, "DOP853 failed: steps={}", result.n_steps);
541 assert_eq!(result.t, t_eval);
542 assert_eq!(result.y[0], y0);
543 }
544
545 #[test]
546 fn test_reports_partial_output_when_t_eval_exceeds_t_span() {
547 let rhs = |_t: f64, y: &[f64], dydt: &mut [f64]| {
548 dydt[0] = -y[0];
549 };
550 let t_eval = [0.0, 0.5, 1.0, 1.5];
551 let result = integrate_dop853(rhs, (0.0, 1.0), &[1.0], &t_eval, &Dop853Options::default());
552 assert!(!result.success, "expected partial output semantics");
553 assert_eq!(result.t, vec![0.0, 0.5, 1.0]);
554 assert_eq!(result.y.len(), result.t.len());
555 }
556}