From d31e1ed3270504f4cb525d02d76ba6a156d8df0d Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Tue, 17 May 2016 16:47:49 +0900 Subject: [PATCH 01/12] rename file --- src/Numerics/Numerics.csproj | 2 +- src/Numerics/OdeSolvers/{OdeSolvers.cs => RungeKutta.cs} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/Numerics/OdeSolvers/{OdeSolvers.cs => RungeKutta.cs} (100%) diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 6369d13e..eb71bc10 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -112,7 +112,7 @@ - + diff --git a/src/Numerics/OdeSolvers/OdeSolvers.cs b/src/Numerics/OdeSolvers/RungeKutta.cs similarity index 100% rename from src/Numerics/OdeSolvers/OdeSolvers.cs rename to src/Numerics/OdeSolvers/RungeKutta.cs From 40e32aaf576c59423f5fc2d4ad63331846688f10 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Tue, 17 May 2016 16:58:48 +0900 Subject: [PATCH 02/12] add Adams-Bashforth First Order(same as Euler method) solver and test --- src/Numerics/Numerics.csproj | 1 + src/Numerics/OdeSolvers/AdamsBashforth.cs | 78 +++++++++++++++++++++++ src/UnitTests/OdeSolvers/OdeSolverTest.cs | 20 ++++++ 3 files changed, 99 insertions(+) create mode 100644 src/Numerics/OdeSolvers/AdamsBashforth.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index eb71bc10..2306ffa1 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -112,6 +112,7 @@ + diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs new file mode 100644 index 00000000..1efee322 --- /dev/null +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -0,0 +1,78 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2016 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// 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) + /// + /// Initial value + /// Start Time + /// End Time + /// Number of subintervals + /// ode model + /// + public static double[] FirstOrder(double y0, double start, double end, int N, Func f) + { + double dt = (end - start) / (N - 1); + double t = start; + double[] y = new double[N]; + y[0] = y0; + for (int i = 1; i < N; i++) + { + y[i] = y0 + dt * f(t, y0); + t += dt; + y0 = y[i]; + } + return y; + } + + public static double[] SecondOrder() + { + return null; + } + + public static double[] ThirdOrder() + { + return null; + } + + public static double[] FourthOrder() + { + return null; + } + } +} diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs index d8b3542b..3d1b18ba 100644 --- a/src/UnitTests/OdeSolvers/OdeSolverTest.cs +++ b/src/UnitTests/OdeSolvers/OdeSolverTest.cs @@ -40,6 +40,26 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers [TestFixture, Category("OdeSolver")] public class OdeSolverTest { + [Test] + public void AB1Test() + { + 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.FirstOrder(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(1, ratio, 0.01);// Check error convergence order + } /// /// Runge-Kutta second order method for first order ODE. /// From 75dd96c9878f0f01078b94c044f9a417ef58449d Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Tue, 17 May 2016 17:10:27 +0900 Subject: [PATCH 03/12] blowing up error when adding second order AB method --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 34 +++++++++++++++++------ src/UnitTests/OdeSolvers/OdeSolverTest.cs | 26 ++++++++++++++++- 2 files changed, 51 insertions(+), 9 deletions(-) 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. /// From 467717436b8c80087f12f5f8d6e10b3123ecf512 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Tue, 17 May 2016 19:22:29 +0900 Subject: [PATCH 04/12] partially fixed (not enough order of convergence) --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 9 ++++++--- src/UnitTests/OdeSolvers/OdeSolverTest.cs | 4 +--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index 825a4b8b..aae22131 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -65,19 +65,22 @@ namespace MathNet.Numerics.OdeSolvers /// Number of subintervals /// ode model /// - public static double[] SecondOrder(double y0, double y1, 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 k1 = f(t, y0); + double k2 = f(t + dt, y0 + dt * k1); + double y1 = y0 + 0.5 * dt * (k1 + k2); y[0] = y0; + t += dt; 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; + y0 = y[i-1]; y1 = y[i]; } return y; diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs index bca1ef38..dfd956fa 100644 --- a/src/UnitTests/OdeSolvers/OdeSolverTest.cs +++ b/src/UnitTests/OdeSolvers/OdeSolverTest.cs @@ -72,9 +72,7 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers 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); + double[] y_t = AdamsBashforth.SecondOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 6)), ode); error = Math.Abs(sol(2) - y_t.Last()); if (oldError != 0) ratio = Math.Log(oldError / error, 2); From 72e9d210406c9cc04dd9814daff65048bc47edf3 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Wed, 18 May 2016 00:16:37 +0900 Subject: [PATCH 05/12] change schem to fill y1 step ( need to convergence test ) --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index aae22131..66ddeb38 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -59,7 +59,6 @@ namespace MathNet.Numerics.OdeSolvers /// Second Order AB Method(Require two initial guesses) /// /// Initial value 1 - /// Initial value 2 /// Start Time /// End Time /// Number of subintervals @@ -70,12 +69,13 @@ namespace MathNet.Numerics.OdeSolvers double dt = (end - start) / (N - 1); double t = start; double[] y = new double[N]; - double k1 = f(t, y0); - double k2 = f(t + dt, y0 + dt * k1); - double y1 = y0 + 0.5 * dt * (k1 + k2); + 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 ); 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)); From 7de883a41a7ca20af27ee2b3b63b03c2fdfdba74 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Wed, 18 May 2016 09:42:53 +0900 Subject: [PATCH 06/12] done Adams-Bashforth second order method to solve ODE --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 14 +++++++------- src/UnitTests/OdeSolvers/OdeSolverTest.cs | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) 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); From 93ddbe6337cd92fbbc787b35c1b51c17d47ed640 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Wed, 18 May 2016 09:54:38 +0900 Subject: [PATCH 07/12] changed initialization method (optimal rk2 -> rk2) --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index 2dbdb5f7..00473028 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -71,9 +71,10 @@ namespace MathNet.Numerics.OdeSolvers double t = start; double[] y = new double[N]; - 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); + double k1 = f(t, y0); + double k2 = f(t + dt, y0 + dt * k1); + double y1 = y0 + 0.5 * dt * (k1 + k2); + y[0] = y0; y[1] = y1; for (int i = 2; i < N; i++) From 98c0a7004b0ef0a05e08fb61470ed5b64ad5b619 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Wed, 18 May 2016 10:25:08 +0900 Subject: [PATCH 08/12] add adams-bashforth 4th order method --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 28 +++++++++++++++++++++-- src/UnitTests/OdeSolvers/OdeSolverTest.cs | 21 +++++++++++++++++ 2 files changed, 47 insertions(+), 2 deletions(-) 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. /// From e1559fe5e1a986d1424f2f760f9f205e786a5af3 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Wed, 18 May 2016 10:50:42 +0900 Subject: [PATCH 09/12] add comment for ab4 --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index 39904dc9..c9b7d2d4 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -92,6 +92,15 @@ namespace MathNet.Numerics.OdeSolvers return null; } + /// + /// Fourth Order AB Method(Require two initial guesses) + /// + /// Initial value 1 + /// Start Time + /// End Time + /// Number of subintervals + /// ode model + /// public static double[] FourthOrder(double y0, double start, double end, int N, Func f) { double dt = (end - start) / (N - 1); From 883bf35c71808d2c912f8950289697ee252046f2 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Wed, 18 May 2016 10:53:45 +0900 Subject: [PATCH 10/12] correct explanation for N --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 6 +++--- src/Numerics/OdeSolvers/RungeKutta.cs | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index c9b7d2d4..42d8ed0d 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -38,7 +38,7 @@ namespace MathNet.Numerics.OdeSolvers /// Initial value /// Start Time /// End Time - /// Number of subintervals + /// Size of output array(the larger, the finer) /// ode model /// public static double[] FirstOrder(double y0, double start, double end, int N, Func f) @@ -62,7 +62,7 @@ namespace MathNet.Numerics.OdeSolvers /// Initial value 1 /// Start Time /// End Time - /// Number of subintervals + /// Size of output array(the larger, the finer) /// ode model /// public static double[] SecondOrder(double y0, double start, double end, int N, Func f) @@ -98,7 +98,7 @@ namespace MathNet.Numerics.OdeSolvers /// Initial value 1 /// Start Time /// End Time - /// Number of subintervals + /// Size of output array(the larger, the finer) /// ode model /// public static double[] FourthOrder(double y0, double start, double end, int N, Func f) diff --git a/src/Numerics/OdeSolvers/RungeKutta.cs b/src/Numerics/OdeSolvers/RungeKutta.cs index dd3ab6fe..bbf82490 100644 --- a/src/Numerics/OdeSolvers/RungeKutta.cs +++ b/src/Numerics/OdeSolvers/RungeKutta.cs @@ -43,7 +43,7 @@ namespace MathNet.Numerics.OdeSolvers /// initial value /// start time /// end time - /// Number of subintervals + /// Size of output array(the larger, the finer) /// ode function /// approximations public static double[] SecondOrder(double y0, double start, double end, int N, Func f) @@ -71,7 +71,7 @@ namespace MathNet.Numerics.OdeSolvers /// initial value /// start time /// end time - /// number of subintervals + /// Size of output array(the larger, the finer) /// ode function /// approximations public static double[] FourthOrder(double y0, double start, double end, int N, Func f) @@ -103,7 +103,7 @@ namespace MathNet.Numerics.OdeSolvers /// initial vector /// start time /// end time - /// number of subintervals + /// Size of output array(the larger, the finer) /// ode function /// approximations public static Vector[] SecondOrder(Vector y0, double start, double end, int N, Func, Vector> f) @@ -130,7 +130,7 @@ namespace MathNet.Numerics.OdeSolvers /// initial vector /// start time /// end time - /// number of subintervals + /// Size of output array(the larger, the finer) /// ode function /// approximations public static Vector[] FourthOrder(Vector y0, double start, double end, int N, Func, Vector> f) From 38d3c39515f43c6567c1323512916c1188ec43fe Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Wed, 18 May 2016 11:04:25 +0900 Subject: [PATCH 11/12] add comments for return object --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Numerics/OdeSolvers/AdamsBashforth.cs b/src/Numerics/OdeSolvers/AdamsBashforth.cs index 42d8ed0d..2519f100 100644 --- a/src/Numerics/OdeSolvers/AdamsBashforth.cs +++ b/src/Numerics/OdeSolvers/AdamsBashforth.cs @@ -40,7 +40,7 @@ namespace MathNet.Numerics.OdeSolvers /// End Time /// Size of output array(the larger, the finer) /// ode model - /// + /// approximation with size N public static double[] FirstOrder(double y0, double start, double end, int N, Func f) { double dt = (end - start) / (N - 1); @@ -64,7 +64,7 @@ namespace MathNet.Numerics.OdeSolvers /// End Time /// Size of output array(the larger, the finer) /// ode model - /// + /// approximation with size N public static double[] SecondOrder(double y0, double start, double end, int N, Func f) { double dt = (end - start) / (N - 1); @@ -100,7 +100,7 @@ namespace MathNet.Numerics.OdeSolvers /// End Time /// Size of output array(the larger, the finer) /// ode model - /// + /// approximation with size N public static double[] FourthOrder(double y0, double start, double end, int N, Func f) { double dt = (end - start) / (N - 1); From 0476e390db141f25cd456143e445a7735cb02d63 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Wed, 18 May 2016 12:52:11 +0900 Subject: [PATCH 12/12] add Adams-Bashforth 3rd order and verify the order --- src/Numerics/OdeSolvers/AdamsBashforth.cs | 41 ++++++++++++++++++++--- src/UnitTests/OdeSolvers/OdeSolverTest.cs | 21 ++++++++++++ 2 files changed, 58 insertions(+), 4 deletions(-) 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() {