diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index 2519f100..0f731537 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -57,7 +57,7 @@ namespace MathNet.Numerics.OdeSolvers } /// - /// Second Order AB Method(Require two initial guesses) + /// Second Order AB Method /// /// Initial value 1 /// Start Time @@ -87,13 +87,46 @@ namespace MathNet.Numerics.OdeSolvers return y; } - public static double[] ThirdOrder() + /// + /// Third Order AB Method + /// + /// Initial value 1 + /// Start Time + /// End Time + /// Size of output array(the larger, the finer) + /// ode model + /// approximation with size N + public static double[] ThirdOrder(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 < 3; 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 = 3; i < N; i++) + { + y[i] = y[i - 1] + dt * (23 * f(t, y[i - 1]) - 16 * f(t - dt, y[i - 2]) + 5 * f(t - 2 * dt, y[i - 3])) / 12.0; + t += dt; + } + return y; } /// - /// Fourth Order AB Method(Require two initial guesses) + /// Fourth Order AB Method /// /// Initial value 1 /// Start Time diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs index 74e3562c..265939a3 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 AB3Test() + { + 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.ThirdOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 8)), 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(3, ratio, 0.01);// Check error convergence order + } + [Test] public void AB4Test() {