diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index 1efee322..825a4b8b 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -27,15 +27,11 @@ // OTHER DEALINGS IN THE SOFTWARE. // using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; namespace MathNet.Numerics.OdeSolvers { public static class AdamsBashforth { - /// /// First Order AB method(same as Forward Euler) /// @@ -59,10 +55,32 @@ namespace MathNet.Numerics.OdeSolvers } return y; } - - public static double[] SecondOrder() + /// + /// Second Order AB Method(Require two initial guesses) + /// + /// Initial value 1 + /// Initial value 2 + /// Start Time + /// End Time + /// Number of subintervals + /// ode model + /// + public static double[] SecondOrder(double y0, double y1, double start, double end, int N, Func f) { - return null; + double dt = (end - start) / (N - 1); + double t = start; + double[] y = new double[N]; + y[0] = y0; + y[1] = y1; + for (int i = 2; i < N; i++) + { + y1 = f(t + dt, y1); + y[i] = y1 + dt * (1.5 * f(t + dt, y1) - 0.5 * f(t, y0)); + t += dt; + y0 = y1; + y1 = y[i]; + } + return y; } public static double[] ThirdOrder() @@ -75,4 +93,4 @@ namespace MathNet.Numerics.OdeSolvers return null; } } -} +} \ No newline at end of file diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs index 3d1b18ba..bca1ef38 100644 --- a/src/UnitTests/OdeSolvers/OdeSolverTest.cs +++ b/src/UnitTests/OdeSolvers/OdeSolverTest.cs @@ -27,9 +27,9 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using MathNet.Numerics.OdeSolvers; using NUnit.Framework; using System; -using MathNet.Numerics.OdeSolvers; using System.Linq; namespace MathNet.Numerics.UnitTests.OdeSolvers @@ -60,6 +60,30 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers } Assert.AreEqual(1, ratio, 0.01);// Check error convergence order } + + [Test] + public void AB2Test() + { + 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; + int N = Convert.ToInt32(Math.Pow(2, k + 6)); + double dt = 2.0 / (N - 1); + double[] y_t = AdamsBashforth.SecondOrder(y0, sol(dt), 0, 2, N, 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(2, ratio, 0.01);// Check error convergence order + } + /// /// Runge-Kutta second order method for first order ODE. ///