diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index 00473028..39904dc9 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/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 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; } } } \ No newline at end of file diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs index 3ac10cd0..74e3562c 100644 --- a/src/UnitTests/OdeSolvers/OdeSolverTest.cs +++ b/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 ode = (t, y) => t + 2 * y * t; + Func 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 + } + /// /// Runge-Kutta second order method for first order ODE. ///