| 19 | return variantToNumberOfStages[variant]; |
| 20 | } |
| 21 | void initializeRungeKuttaScheme(RungeKuttaVariant variant, |
| 22 | int& numberOfStages, |
| 23 | Eigen::MatrixXd& a, |
| 24 | Eigen::VectorXd& b, |
| 25 | Eigen::VectorXd& c) { |
| 26 | numberOfStages = getNumberOfStages(variant); |
| 27 | |
| 28 | // Initialize coefficients |
| 29 | a = Eigen::MatrixXd(numberOfStages, numberOfStages); |
| 30 | b = Eigen::VectorXd(numberOfStages); |
| 31 | c = Eigen::VectorXd(numberOfStages); |
| 32 | |
| 33 | a.setZero(); |
| 34 | b.setZero(); |
| 35 | c.setZero(); |
| 36 | |
| 37 | switch (variant) { |
| 38 | case RungeKuttaVariant::RK4: |
| 39 | // The classical RK4 |
| 40 | a(1, 0) = 1.0 / 2.0; |
| 41 | a(2, 0) = 0.0; |
| 42 | a(2, 1) = 1.0 / 2.0; |
| 43 | a(3, 0) = 0.0; |
| 44 | a(3, 1) = 0.0; |
| 45 | a(3, 2) = 1.0; |
| 46 | |
| 47 | b(0) = 1.0 / 6.0; |
| 48 | b(1) = 1.0 / 3.0; |
| 49 | b(2) = 1.0 / 3.0; |
| 50 | b(3) = 1.0 / 6.0; |
| 51 | |
| 52 | c(0) = 0.0; |
| 53 | c(1) = 1.0 / 2.0; |
| 54 | c(2) = 1.0 / 2.0; |
| 55 | c(3) = 1.0; |
| 56 | break; |
| 57 | case RungeKuttaVariant::RK438: |
| 58 | // The also classical 3/8 rule |
| 59 | a(1, 0) = 1.0 / 3.0; |
| 60 | a(2, 0) = -1.0 / 3.0; |
| 61 | a(2, 1) = 1.0; |
| 62 | a(3, 0) = 1.0; |
| 63 | a(3, 1) = -1.0; |
| 64 | a(3, 2) = 1.0; |
| 65 | |
| 66 | b(0) = 1.0 / 8.0; |
| 67 | b(1) = 3.0 / 8.0; |
| 68 | b(2) = 3.0 / 8.0; |
| 69 | b(3) = 1.0 / 8.0; |
| 70 | |
| 71 | c(0) = 0.0; |
| 72 | c(1) = 1.0 / 3.0; |
| 73 | c(2) = 2.0 / 3.0; |
| 74 | c(3) = 1.0; |
| 75 | break; |
| 76 | case RungeKuttaVariant::RK4Ralston: |
| 77 | // Ralston's RK4, minimized truncation error. Coeffs stolen from: |
| 78 | // https://github.com/SciML/DiffEqDevTools.jl/blob/b5aca9330cd1a1b6ffbdbdf33a7ea037f7b53699/src/ode_tableaus.jl#L235 |
no test coverage detected