Add Dormand-Prince 8(5)3 method

No interpolant yet though.

Signed-off-by: Lamont Granquist <lamont@scriptkiddie.org>

MechJebLib/MechJebLib.csproj

@@ -78,6 +78,7 @@
         <Compile Include="ODE\AbstractRungeKutta.cs"/>
         <Compile Include="ODE\BS3.cs"/>
         <Compile Include="ODE\DP5.cs"/>
+        <Compile Include="ODE\DP8.cs"/>
         <Compile Include="ODE\Event.cs"/>
         <Compile Include="Primitives\Dual.cs"/>
         <Compile Include="Primitives\DualV3.cs"/>

MechJebLib/ODE/AbstractRungeKutta.cs

@@ -5,7 +5,6 @@
 
 using System;
 using System.Collections.Generic;
-using JetBrains.Annotations;
 using MechJebLib.Primitives;
 using static MechJebLib.Utils.Statics;
 using static System.Math;
@@ -39,8 +38,6 @@ public double Beta
 
         protected override (double, double) Step(IVPFunc f)
         {
-            using var err = Vn.Rent(N);
-
             bool previouslyRejected = false;
 
             while (true)
@@ -52,17 +49,13 @@ protected override (double, double) Step(IVPFunc f)
                 else if (Habs < MinStep)
                     Habs = MinStep;
 
-                RKStep(f, err);
+                RKStep(f);
 
-                double errorNorm = ScaledErrorNorm(err);
+                double errorNorm = ScaledErrorNorm();
 
                 if (errorNorm < 1)
                 {
-                    double factor;
-                    if (errorNorm == 0)
-                        factor = MAX_FACTOR;
-                    else
-                        factor = Min(MAX_FACTOR, SAFETY * Pow(errorNorm, -_alpha) * Pow(_lastErrorNorm, Beta));
+                    double factor = errorNorm == 0 ? MAX_FACTOR : Min(MAX_FACTOR, SAFETY * Pow(errorNorm, -_alpha) * Pow(_lastErrorNorm, Beta));
 
                     if (previouslyRejected)
                         factor = Min(1.0, factor);
@@ -129,22 +122,8 @@ protected override void Cleanup()
             K.Clear();
         }
 
-        protected abstract void RKStep(IVPFunc f, Vn err);
-
-        [UsedImplicitly]
-        protected virtual double ScaledErrorNorm(Vn err)
-        {
-            int n = err.Count;
-
-            double error = 0.0;
+        protected abstract void RKStep(IVPFunc f);
 
-            for (int i = 0; i < n; i++)
-            {
-                double scale = Atol + Rtol * Max(Abs(Y[i]), Abs(Ynew[i]));
-                error += Powi(err[i] / scale, 2);
-            }
-
-            return Sqrt(error / n);
-        }
+        protected abstract double ScaledErrorNorm();
     }
 }

MechJebLib/ODE/BS3.cs

@@ -7,6 +7,8 @@
 using System.Collections.Generic;
 using MechJebLib.Functions;
 using MechJebLib.Primitives;
+using static MechJebLib.Utils.Statics;
+using static System.Math;
 
 namespace MechJebLib.ODE
 {
@@ -40,7 +42,7 @@ public class BS3 : AbstractRungeKutta
 
         #endregion
 
-        protected override void RKStep(IVPFunc f, Vn err)
+        protected override void RKStep(IVPFunc f)
         {
             double h = Habs * Direction;
 
@@ -58,10 +60,27 @@ protected override void RKStep(IVPFunc f, Vn err)
                 Ynew[i] = Y[i] + h * (A41 * Dy[i] + A42 * K[2][i] + A43 * K[3][i]);
             f(Ynew, T + h, K[4]);
 
+
+
+            K[4].CopyTo(Dynew);
+        }
+
+        protected override double ScaledErrorNorm()
+        {
+            using var err = Vn.Rent(N);
+
             for (int i = 0; i < N; i++)
                 err[i] = Dy[i] * E1 + K[2][i] * E2 + K[3][i] * E3 + K[4][i] * E4;
 
-            K[4].CopyTo(Dynew);
+            double error = 0.0;
+
+            for (int i = 0; i < N; i++)
+            {
+                double scale = Atol + Rtol * Max(Abs(Y[i]), Abs(Ynew[i]));
+                error += Powi(err[i] / scale, 2);
+            }
+
+            return Sqrt(error / N);
         }
 
         protected override void InitInterpolant()

MechJebLib/ODE/DP5.cs

@@ -6,6 +6,8 @@
 using System;
 using System.Collections.Generic;
 using MechJebLib.Primitives;
+using static MechJebLib.Utils.Statics;
+using static System.Math;
 
 // ReSharper disable CompareOfFloatsByEqualityOperator
 namespace MechJebLib.ODE
@@ -94,7 +96,7 @@ public class DP5 : AbstractRungeKutta
 
         #endregion
 
-        protected override void RKStep(IVPFunc f, Vn err)
+        protected override void RKStep(IVPFunc f)
         {
             double h = Habs * Direction;
 
@@ -125,10 +127,26 @@ protected override void RKStep(IVPFunc f, Vn err)
 
             f(Ynew, T + h, K[7]);
 
+
+            K[7].CopyTo(Dynew);
+        }
+
+        protected override double ScaledErrorNorm()
+        {
+            using var err = Vn.Rent(N);
+
             for (int i = 0; i < N; i++)
                 err[i] = K[1][i] * E1 + K[3][i] * E3 + K[4][i] * E4 + K[5][i] * E5 + K[6][i] * E6 + K[7][i] * E7;
 
-            K[7].CopyTo(Dynew);
+            double error = 0.0;
+
+            for (int i = 0; i < N; i++)
+            {
+                double scale = Atol + Rtol * Max(Abs(Y[i]), Abs(Ynew[i]));
+                error += Powi(err[i] / scale, 2);
+            }
+
+            return Sqrt(error / N);
         }
 
         protected override void InitInterpolant()
@@ -153,9 +171,7 @@ protected override void Interpolate(double x, Vn yout)
             double bs7 = (1.0 - s) * s * (B71 + B72 * s + B73 * s2);
 
             for (int i = 0; i < N; i++)
-            {
                 yout[i] = Y[i] + h * s * (bs1 * K[1][i] + bs3 * K[3][i] + bs4 * K[4][i] + bs5 * K[5][i] + bs6 * K[6][i] + bs7 * K[7][i]);
-            }
         }
     }
 }