From 1bbe6c897aa69622288d430504ec45c6c8169c1b Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Fri, 13 May 2016 20:18:46 +0900 Subject: [PATCH 1/6] add RK2 and its test to verify RK2 erro convergence order --- src/Numerics/Numerics.csproj | 1 + src/Numerics/OdeSolvers/OdeSolvers.cs | 59 ++++++++++++++++++++++ src/UnitTests/OdeSolvers/OdeSolverTest.cs | 61 +++++++++++++++++++++++ src/UnitTests/UnitTests.csproj | 1 + 4 files changed, 122 insertions(+) create mode 100644 src/Numerics/OdeSolvers/OdeSolvers.cs create mode 100644 src/UnitTests/OdeSolvers/OdeSolverTest.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index b6004e90..6369d13e 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -112,6 +112,7 @@ + diff --git a/src/Numerics/OdeSolvers/OdeSolvers.cs b/src/Numerics/OdeSolvers/OdeSolvers.cs new file mode 100644 index 00000000..5bab7436 --- /dev/null +++ b/src/Numerics/OdeSolvers/OdeSolvers.cs @@ -0,0 +1,59 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-2013 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 MathNet.Numerics.Properties; + +namespace MathNet.Numerics.OdeSolvers +{ + /// + /// ODE Solver Algorithms + /// + public static class RungeKutta + { + public static double[] SecondOrder(double y0, double start, double end, int N, Func f) + { + double dt = (end - start) / (N - 1); + double k1 = 0; + double k2 = 0; + double t = start; + double[] y = new double[N]; + y[0] = y0; + for (int i = 1; i < N; i++) + { + k1 = f(t, y0); + k2 = f(t + dt, y0 + k1 * dt); + y[i] = y0 + dt * 0.5 * (k1 + k2); + t += dt; + y0 = y[i]; + } + return y; + } + } +} \ No newline at end of file diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs new file mode 100644 index 00000000..342de719 --- /dev/null +++ b/src/UnitTests/OdeSolvers/OdeSolverTest.cs @@ -0,0 +1,61 @@ +// +// 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. +// + +extern alias NUnitFramework; + +using System; +using MathNet.Numerics.OdeSolvers; +using NUnitFramework.NUnit.Framework; +using System.Linq; + +namespace MathNet.Numerics.UnitTests.OdeSolvers +{ + /// + /// ODE Solver tests. + /// + [TestFixture, Category("OdeSolver")] + public class OdeSolverTest + { + /// + /// Runge-Kutta second order method for first order ODE. + /// + [Test] + public void RK2Test() + { + Func ode = (t, y) => t + 2 * y * t; + Func sol = (t) => 0.5 * (Math.Exp(t * t) - 1); + for (int k = 0; k < 4; k++) + { + double y0 = 0; + double[] y_t = RungeKutta.SecondOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 6)), ode); + Console.WriteLine(Math.Abs(sol(2) - y_t.Last())); + } + } + } +} \ No newline at end of file diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 735d4a80..e07050ea 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -362,6 +362,7 @@ + From 02a0af61a236ebf68a47e6b15d0d33945a7a2bd7 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Fri, 13 May 2016 20:23:46 +0900 Subject: [PATCH 2/6] add RK4 and its test to verify convergence order --- src/Numerics/OdeSolvers/OdeSolvers.cs | 23 +++++++++++++++++++++++ src/UnitTests/OdeSolvers/OdeSolverTest.cs | 16 ++++++++++++++++ 2 files changed, 39 insertions(+) diff --git a/src/Numerics/OdeSolvers/OdeSolvers.cs b/src/Numerics/OdeSolvers/OdeSolvers.cs index 5bab7436..c630710c 100644 --- a/src/Numerics/OdeSolvers/OdeSolvers.cs +++ b/src/Numerics/OdeSolvers/OdeSolvers.cs @@ -55,5 +55,28 @@ namespace MathNet.Numerics.OdeSolvers } return y; } + + public static double[] FourthOrder(double y0, double start, double end, int N, Func f) + { + double dt = (end - start) / (N - 1); + double k1 = 0; + double k2 = 0; + double k3 = 0; + double k4 = 0; + double t = start; + double[] y = new double[N]; + y[0] = y0; + for (int i = 1; i < N; i++) + { + k1 = f(t, y0); + k2 = f(t + dt / 2, y0 + k1 * dt / 2); + k3 = f(t + dt / 2, y0 + k2 * dt / 2); + k4 = f(t + dt, y0 + k3 * dt); + y[i] = y0 + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4); + t += dt; + y0 = y[i]; + } + return y; + } } } \ No newline at end of file diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs index 342de719..f57b855c 100644 --- a/src/UnitTests/OdeSolvers/OdeSolverTest.cs +++ b/src/UnitTests/OdeSolvers/OdeSolverTest.cs @@ -57,5 +57,21 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers Console.WriteLine(Math.Abs(sol(2) - y_t.Last())); } } + + /// + /// Runge-Kutta second order method for first order ODE. + /// + [Test] + public void RK4Test() + { + Func ode = (t, y) => t + 2 * y * t; + Func sol = (t) => 0.5 * (Math.Exp(t * t) - 1); + for (int k = 0; k < 4; k++) + { + double y0 = 0; + double[] y_t = RungeKutta.FourthOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 6)), ode); + Console.WriteLine(Math.Abs(sol(2) - y_t.Last())); + } + } } } \ No newline at end of file From b204ac0dcf97d7cb90975318232ae1861edc701c Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Fri, 13 May 2016 20:25:15 +0900 Subject: [PATCH 3/6] add RK2 and RK4 to solve ODE System --- src/Numerics/OdeSolvers/OdeSolvers.cs | 39 +++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/src/Numerics/OdeSolvers/OdeSolvers.cs b/src/Numerics/OdeSolvers/OdeSolvers.cs index c630710c..0cbb877a 100644 --- a/src/Numerics/OdeSolvers/OdeSolvers.cs +++ b/src/Numerics/OdeSolvers/OdeSolvers.cs @@ -29,6 +29,7 @@ using System; using MathNet.Numerics.Properties; +using MathNet.Numerics.LinearAlgebra; namespace MathNet.Numerics.OdeSolvers { @@ -78,5 +79,43 @@ namespace MathNet.Numerics.OdeSolvers } return y; } + + public static Vector[] SecondOrder(Vector y0, double start, double end, int N, Func, Vector> f) + { + double dt = (end - start) / (N - 1); + Vector k1, k2; + Vector[] y = new Vector[N]; + double t = start; + y[0] = y0; + for (int i = 1; i < N; i++) + { + k1 = f(t, y0); + k2 = f(t, y0 + k1 + dt); + y[i] = y0 + 0.5 * (k1 + k2); + t += dt; + y0 = y[i]; + } + return y; + } + + public static Vector[] FourthOrder(Vector y0, double start, double end, int N, Func, Vector> f) + { + double dt = (end - start) / (N - 1); + Vector k1, k2, k3, k4; + Vector[] y = new Vector[N]; + double t = start; + y[0] = y0; + for (int i = 1; i < N; i++) + { + k1 = f(t, y0); + k2 = f(t + dt / 2, y0 + k1 * dt / 2); + k3 = f(t + dt / 2, y0 + k2 * dt / 2); + k4 = f(t + dt, y0 + k3 * dt); + y[i] = y0 + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4); + t += dt; + y0 = y[i]; + } + return y; + } } } \ No newline at end of file From 72ed2931fe775b4d0bfe7afbbe7f90e27d3de7d0 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Fri, 13 May 2016 20:31:52 +0900 Subject: [PATCH 4/6] typo --- src/Numerics/OdeSolvers/OdeSolvers.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Numerics/OdeSolvers/OdeSolvers.cs b/src/Numerics/OdeSolvers/OdeSolvers.cs index 0cbb877a..136f7b7b 100644 --- a/src/Numerics/OdeSolvers/OdeSolvers.cs +++ b/src/Numerics/OdeSolvers/OdeSolvers.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2013 Math.NET +// 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 From b2c3fcaf3272e6e56d2e3a3b643d0fffcc99d446 Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Mon, 16 May 2016 19:19:22 +0900 Subject: [PATCH 5/6] add comments and check the error convergence --- src/Numerics/OdeSolvers/OdeSolvers.cs | 39 +++++++++++++++++++++-- src/UnitTests/OdeSolvers/OdeSolverTest.cs | 22 +++++++++++-- 2 files changed, 56 insertions(+), 5 deletions(-) diff --git a/src/Numerics/OdeSolvers/OdeSolvers.cs b/src/Numerics/OdeSolvers/OdeSolvers.cs index 136f7b7b..dd3ab6fe 100644 --- a/src/Numerics/OdeSolvers/OdeSolvers.cs +++ b/src/Numerics/OdeSolvers/OdeSolvers.cs @@ -27,9 +27,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; -using MathNet.Numerics.Properties; using MathNet.Numerics.LinearAlgebra; +using System; namespace MathNet.Numerics.OdeSolvers { @@ -38,6 +37,15 @@ namespace MathNet.Numerics.OdeSolvers /// public static class RungeKutta { + /// + /// Second Order Runge-Kutta method + /// + /// initial value + /// start time + /// end time + /// Number of subintervals + /// ode function + /// approximations public static double[] SecondOrder(double y0, double start, double end, int N, Func f) { double dt = (end - start) / (N - 1); @@ -57,6 +65,15 @@ namespace MathNet.Numerics.OdeSolvers return y; } + /// + /// Fourth Order Runge-Kutta method + /// + /// initial value + /// start time + /// end time + /// number of subintervals + /// ode function + /// approximations public static double[] FourthOrder(double y0, double start, double end, int N, Func f) { double dt = (end - start) / (N - 1); @@ -80,6 +97,15 @@ namespace MathNet.Numerics.OdeSolvers return y; } + /// + /// Second Order Runge-Kutta to solve ODE SYSTEM + /// + /// initial vector + /// start time + /// end time + /// number of subintervals + /// ode function + /// approximations public static Vector[] SecondOrder(Vector y0, double start, double end, int N, Func, Vector> f) { double dt = (end - start) / (N - 1); @@ -98,6 +124,15 @@ namespace MathNet.Numerics.OdeSolvers return y; } + /// + /// Fourth Order Runge-Kutta to solve ODE SYSTEM + /// + /// initial vector + /// start time + /// end time + /// number of subintervals + /// ode function + /// approximations public static Vector[] FourthOrder(Vector y0, double start, double end, int N, Func, Vector> f) { double dt = (end - start) / (N - 1); diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs index f57b855c..7fb828b9 100644 --- a/src/UnitTests/OdeSolvers/OdeSolverTest.cs +++ b/src/UnitTests/OdeSolvers/OdeSolverTest.cs @@ -50,28 +50,44 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers { 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 = RungeKutta.SecondOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 6)), ode); - Console.WriteLine(Math.Abs(sol(2) - y_t.Last())); + 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. + /// Runge-Kutta fourth order method for first order ODE. /// [Test] public void RK4Test() { 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 = RungeKutta.FourthOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 6)), ode); - Console.WriteLine(Math.Abs(sol(2) - y_t.Last())); + 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 } } } \ No newline at end of file From dc75e1589c97a2e013897708d6e8999687a0b8ca Mon Sep 17 00:00:00 2001 From: Yoonku Hwang Date: Mon, 16 May 2016 20:15:32 +0900 Subject: [PATCH 6/6] change NUnit using lines --- src/UnitTests/OdeSolvers/OdeSolverTest.cs | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/UnitTests/OdeSolvers/OdeSolverTest.cs b/src/UnitTests/OdeSolvers/OdeSolverTest.cs index 7fb828b9..d8b3542b 100644 --- a/src/UnitTests/OdeSolvers/OdeSolverTest.cs +++ b/src/UnitTests/OdeSolvers/OdeSolverTest.cs @@ -27,11 +27,9 @@ // OTHER DEALINGS IN THE SOFTWARE. // -extern alias NUnitFramework; - +using NUnit.Framework; using System; using MathNet.Numerics.OdeSolvers; -using NUnitFramework.NUnit.Framework; using System.Linq; namespace MathNet.Numerics.UnitTests.OdeSolvers