RungeKutta/src/RKOneCashKarpStep.icc

0001 template <typename T, int N>
0002 std::pair<typename RKOneCashKarpStep<T, N>::Vector, T> RKOneCashKarpStep<T, N>::operator()(
0003     Scalar x, const Vector& v, const RKDerivative<T, N>& deriv, const RKDistance<T, N>& dist, Scalar step) {
0004   const Scalar a2 = 0.2, a3 = 0.3, a4 = 0.6, a5 = 1., a6 = 7. / 8.;
0005   const Scalar b21 = 0.2;
0006   const Scalar b31 = 3. / 40., b32 = 9. / 40.;
0007   const Scalar b41 = 0.3, b42 = -0.9, b43 = 1.2;
0008   const Scalar b51 = -11. / 54., b52 = 5. / 2., b53 = -70. / 27., b54 = 35. / 27.;
0009   const Scalar b61 = 1631. / 55296., b62 = 175. / 512., b63 = 575. / 13824., b64 = 44275. / 110592., b65 = 253. / 4096.;
0010   const Scalar c1 = 37. / 378., c3 = 250. / 621., c4 = 125. / 594., c6 = 512. / 1771.;
0011   // removed unused variables c2=0, c5=0 MM 25/6/07
0012   const Scalar d1 = 2825. / 27648., d3 = 18575. / 48384., d4 = 13525. / 55296., d5 = 277. / 14336., d6 = 0.25;
0013   // reomved unused variable d2=0
0014
0015   Vector k1 = step * deriv(x, v);
0016   Vector k2 = step * deriv(x + a2 * step, v + b21 * k1);
0017   Vector k3 = step * deriv(x + a3 * step, v + b31 * k1 + b32 * k2);
0018   Vector k4 = step * deriv(x + a4 * step, v + b41 * k1 + b42 * k2 + b43 * k3);
0019   Vector k5 = step * deriv(x + a5 * step, v + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4);
0020   Vector k6 = step * deriv(x + a6 * step, v + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5);
0021
0022   Vector r5 = v + c1 * k1 + c3 * k3 + c4 * k4 + c6 * k6;
0023   Vector r4 = v + d1 * k1 + d3 * k3 + d4 * k4 + d5 * k5 + d6 * k6;
0024   return std::pair<Vector, T>(r5, dist(r4, r5, x + step));
0025 }
0026
0027 /* array implementation may be faster, but is tricky (needs triangular matrix etc.)
0028 {
0029   Scalar a[6] = {0, 0.2, 0.3, 0.6, 1., 7./8.};
0030
0031   k[1] = step*deriv( x, v);
0032   for (int i=1; i<=6; i++) {
0033     Vector arg = v;
0034     for (int j=1; j<i; j++) v.increment( k[j], b[i][j]); // v += b[i][j]*k[j];
0035     k[i] = step*deriv( x+a[i], arg);
0036   }
0037
0038 }
0039 */