diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index 66ddeb38..2dbdb5f7 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -55,6 +55,7 @@ namespace MathNet.Numerics.OdeSolvers } return y; } + /// /// Second Order AB Method(Require two initial guesses) /// @@ -64,23 +65,22 @@ namespace MathNet.Numerics.OdeSolvers /// Number of subintervals /// ode model /// - public static double[] SecondOrder(double y0, double start, double end, int N, Func f) + public static double[] SecondOrder(double y0, double start, double end, int N, Func f) { double dt = (end - start) / (N - 1); double t = start; double[] y = new double[N]; - double fold = f(t, y0); - double ytilde = y0 + (2*dt/3) * fold; - double y1 = y0 + (dt/4) * fold + (3*dt/4) * f(t + (2*dt/3), ytilde ); + + double fold = f(t, y0); + double ytilde = y0 + 2.0 * dt / 3.0 * fold; + double y1 = y0 + dt / 4.0 * fold + 3.0 * dt / 4 * f(t + 2.0 / 3.0 * dt, ytilde); y[0] = y0; - t += dt; y[1] = y1; - for (int i = 2; i < N; i++) { y[i] = y1 + dt * (1.5 * f(t + dt, y1) - 0.5 * f(t, y0)); t += dt; - y0 = y[i-1]; + y0 = y[i - 1]; y1 = y[i]; } return y; diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs index dfd956fa..3ac10cd0 100644 --- a/src/UnitTests/OdeSolvers/OdeSolverTest.cs +++ b/src/UnitTests/OdeSolvers/OdeSolverTest.cs @@ -72,7 +72,7 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers for (int k = 0; k < 4; k++) { double y0 = 0; - double[] y_t = AdamsBashforth.SecondOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 6)), ode); + double[] y_t = AdamsBashforth.SecondOrder(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);