add adams-bashforth 4th order method

src/Numerics/OdeSolvers/AdamsBashforth.cs

@@ -92,9 +92,33 @@ namespace MathNet.Numerics.OdeSolvers
             return null;
         }
 
-        public static double[] FourthOrder()
+        public static double[] FourthOrder(double y0, double start, double end, int N, Func<double, double, double> f)
         {
-            return null;
+            double dt = (end - start) / (N - 1);
+            double t = start;
+            double[] y = new double[N];
+
+            double k1 = 0;
+            double k2 = 0;
+            double k3 = 0;
+            double k4 = 0;
+            y[0] = y0;
+            for (int i = 1; i < 4; i++)
+            {
+                k1 = dt * f(t, y0);
+                k2 = dt * f(t + dt / 2, y0 + k1 / 2);
+                k3 = dt * f(t + dt / 2, y0 + k2 / 2);
+                k4 = dt * f(t + dt, y0 + k3);
+                y[i] = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
+                t += dt;
+                y0 = y[i];
+            }
+            for (int i = 4; i < N; i++)
+            {
+                y[i] = y[i - 1] + dt * (55 * f(t, y[i - 1]) - 59 * f(t - dt, y[i - 2]) + 37 * f(t - 2 * dt, y[i - 3]) - 9 * f(t - 3 * dt, y[i - 4])) / 24.0;
+                t += dt;
+            }
+            return y;
         }
     }
 }

src/UnitTests/OdeSolvers/OdeSolverTest.cs

@@ -82,6 +82,27 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers
             Assert.AreEqual(2, ratio, 0.01);// Check error convergence order
         }
 
+        [Test]
+        public void AB4Test()
+        {
+            Func<double, double, double> ode = (t, y) => t + 2 * y * t;
+            Func<double, double> sol = (t) => 0.5 * (Math.Exp(t * t) - 1);
+            double ratio = double.NaN;
+            double error = 0;
+            double oldError = 0;
+            for (int k = 0; k < 4; k++)
+            {
+                double y0 = 0;
+                double[] y_t = AdamsBashforth.FourthOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 9)) + 1, ode);
+                error = Math.Abs(sol(2) - y_t.Last());
+                if (oldError != 0)
+                    ratio = Math.Log(oldError / error, 2);
+                oldError = error;
+                Console.WriteLine(string.Format("{0}, {1}", error, ratio));
+            }
+            Assert.AreEqual(4, ratio, 0.01);// Check error convergence order
+        }
+
         /// <summary>
         /// Runge-Kutta second order method for first order ODE.
         /// </summary>