numeris 0.5.2

Pure-Rust numerical algorithms library — high performance with SIMD support while also supporting no-std for embedded and WASM targets.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100

//! Verner 9(8) Efficient (no interpolation) — 16 stages, order 9(8).
//!
//! Same core 16 stages as [`RKV98`] but without the 5 extra interpolation
//! stages. Use this when you only need the final answer.
//!
//! Source: <https://www.sfu.ca/~jverner/RKV98.IIa.Efficient.000000349.081209.FLOAT6040OnWeb>

use super::adaptive::RKAdaptive;
#[cfg(feature = "std")]
use super::{OdeError, Solution};
#[cfg(feature = "std")]
use crate::traits::FloatScalar;
#[cfg(feature = "std")]
use crate::Matrix;

/// Verner 9(8) Efficient (no interpolation) — 16 stages, order 9(8).
pub struct RKV98NoInterp;

const BHAT: [f64; 16] = [
    0.019_969_965_148_867_73,
    0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    2.191_499_304_949_33,
    0.088_570_718_482_084_38,
    0.114_056_023_486_596_56,
    0.253_316_380_534_510_7,
    -2.056_564_386_240_941,
    0.340_809_679_901_312,
    0.0,
    0.0,
    0.048_342_313_738_239_585,
];

impl RKAdaptive<16, 1> for RKV98NoInterp {
    const ORDER: usize = 9;
    const FSAL: bool = false;

    const C: [f64; 16] = [
        0.0, 0.034_62, 0.097_024_350_638_780_44, 0.145_536_525_958_170_68,
        0.561, 0.229_007_911_590_485, 0.544_992_088_409_515, 0.645,
        0.483_75, 0.067_57, 0.25, 0.659_065_061_873_099_9,
        0.820_6, 0.901_2, 1.0, 1.0,
    ];

    const B: [f64; 16] = [
        0.014_611_976_858_423_152,
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
        -0.391_521_186_233_133_9,
        0.231_093_250_028_950_65,
        0.127_476_676_999_285_25,
        0.224_643_417_620_415_8,
        0.568_435_268_974_851_3,
        0.058_258_715_572_158_275,
        0.136_431_740_348_221_56,
        0.030_570_139_830_827_976,
        0.0,
    ];

    const BERR: [f64; 16] = {
        let mut berr = [0.0; 16];
        let mut ix: usize = 0;
        while ix < 16 {
            berr[ix] = RKV98NoInterp::B[ix] - BHAT[ix];
            ix += 1;
        }
        berr
    };

    const BI: [[f64; 1]; 16] = [
        [1.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],
        [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0],
    ];

    #[rustfmt::skip]
    const A: [[f64; 16]; 16] = [
        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [0.034_62, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [-0.038_933_543_885_728_734, 0.135_957_894_524_509_18, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [0.036_384_131_489_542_67, 0.0, 0.109_152_394_468_628, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [2.025_763_914_393_97, 0.0, -7.638_023_836_496_292, 6.173_259_922_102_322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [0.051_122_755_894_060_61, 0.0, 0.0, 0.177_082_379_455_502_15, 0.000_802_776_240_922_250_2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [0.131_600_635_797_521_63, 0.0, 0.0, -0.295_727_625_266_963_67, 0.087_813_780_356_429_52, 0.621_305_297_522_527_5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [0.071_666_666_666_666_67, 0.0, 0.0, 0.0, 0.0, 0.330_553_357_891_531_95, 0.242_779_975_441_801_38, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [0.071_806_640_625, 0.0, 0.0, 0.0, 0.0, 0.329_438_028_322_817_7, 0.116_519_002_927_182_29, -0.034_013_671_875, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [0.048_367_576_463_406_47, 0.0, 0.0, 0.0, 0.0, 0.039_289_899_256_761_64, 0.105_474_094_589_034_46, -0.021_438_652_846_483_126, -0.104_122_917_462_719_44, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [-0.026_645_614_872_014_785, 0.0, 0.0, 0.0, 0.0, 0.033_333_333_333_333_33, -0.163_107_224_487_246_7, 0.033_960_816_841_277_615, 0.157_231_941_381_462_6, 0.215_226_747_803_187_96, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
        [0.036_890_092_487_086_225, 0.0, 0.0, 0.0, 0.0, -0.146_518_157_672_554_3, 0.224_257_776_817_202_44, 0.022_944_057_170_660_725, -0.003_585_005_290_572_876, 0.086_692_233_164_443_85, 0.438_384_065_196_833_76, 0.0, 0.0, 0.0, 0.0, 0.0],
        [-0.486_601_221_511_334_06, 0.0, 0.0, 0.0, 0.0, -6.304_602_650_282_853, -0.281_245_618_289_472_6, -2.679_019_236_219_849_3, 0.518_815_663_924_157_6, 1.365_353_187_603_341_8, 5.885_091_088_503_946_5, 2.802_808_786_272_062_8, 0.0, 0.0, 0.0, 0.0],
        [0.418_536_745_775_347_16, 0.0, 0.0, 0.0, 0.0, 6.724_547_581_906_459, -0.425_444_280_164_611_8, 3.343_279_153_001_265_8, 0.617_081_663_117_537_8, -0.929_966_123_939_932_8, -6.099_948_804_751_011, -3.002_206_187_889_399, 0.255_320_252_944_344_6, 0.0, 0.0, 0.0],
        [-0.779_374_086_122_884_6, 0.0, 0.0, 0.0, 0.0, -13.937_342_538_107_776, 1.252_048_853_379_357_2, -14.691_500_408_016_87, -0.494_705_058_533_141_7, 2.242_974_909_146_236_8, 13.367_893_803_828_643, 14.396_650_486_650_687, -0.797_581_333_177_68, 0.440_935_370_953_427_8, 0.0, 0.0],
        [2.058_051_337_466_886_3, 0.0, 0.0, 0.0, 0.0, 22.357_937_727_968_032, 0.909_498_109_975_563_4, 35.891_100_982_402_64, -3.442_515_027_624_453_6, -4.865_481_358_036_368_5, -18.909_803_813_543_427, -34.263_544_480_304_52, 1.264_756_521_695_642_7, 0.0, 0.0, 0.0],
    ];

    #[cfg(feature = "std")]
    fn interpolate<T: FloatScalar, const M: usize, const N: usize>(
        _t_interp: T,
        _sol: &Solution<T, M, N>,
    ) -> Result<Matrix<T, M, N>, OdeError> {
        Err(OdeError::InterpNotImplemented)
    }
}