From b2c3fcaf3272e6e56d2e3a3b643d0fffcc99d446 Mon Sep 17 00:00:00 2001
From: Yoonku Hwang <yz0624@gmail.com>
Date: Mon, 16 May 2016 19:19:22 +0900
Subject: [PATCH] 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.
 // </copyright>
 
-using System;
-using MathNet.Numerics.Properties;
 using MathNet.Numerics.LinearAlgebra;
+using System;
 
 namespace MathNet.Numerics.OdeSolvers
 {
@@ -38,6 +37,15 @@ namespace MathNet.Numerics.OdeSolvers
     /// </summary>
     public static class RungeKutta
     {
+        /// <summary>
+        /// Second Order Runge-Kutta method
+        /// </summary>
+        /// <param name="y0">initial value</param>
+        /// <param name="start">start time</param>
+        /// <param name="end">end time</param>
+        /// <param name="N">Number of subintervals</param>
+        /// <param name="f">ode function</param>
+        /// <returns>approximations</returns>
         public static double[] SecondOrder(double y0, double start, double end, int N, Func<double, double, double> f)
         {
             double dt = (end - start) / (N - 1);
@@ -57,6 +65,15 @@ namespace MathNet.Numerics.OdeSolvers
             return y;
         }
 
+        /// <summary>
+        /// Fourth Order Runge-Kutta method
+        /// </summary>
+        /// <param name="y0">initial value</param>
+        /// <param name="start">start time</param>
+        /// <param name="end">end time</param>
+        /// <param name="N">number of subintervals</param>
+        /// <param name="f">ode function</param>
+        /// <returns>approximations</returns>
         public static double[] FourthOrder(double y0, double start, double end, int N, Func<double, double, double> f)
         {
             double dt = (end - start) / (N - 1);
@@ -80,6 +97,15 @@ namespace MathNet.Numerics.OdeSolvers
             return y;
         }
 
+        /// <summary>
+        /// Second Order Runge-Kutta to solve ODE SYSTEM
+        /// </summary>
+        /// <param name="y0">initial vector</param>
+        /// <param name="start">start time</param>
+        /// <param name="end">end time</param>
+        /// <param name="N">number of subintervals</param>
+        /// <param name="f">ode function</param>
+        /// <returns>approximations</returns>
         public static Vector<double>[] SecondOrder(Vector<double> y0, double start, double end, int N, Func<double, Vector<double>, Vector<double>> f)
         {
             double dt = (end - start) / (N - 1);
@@ -98,6 +124,15 @@ namespace MathNet.Numerics.OdeSolvers
             return y;
         }
 
+        /// <summary>
+        /// Fourth Order Runge-Kutta to solve ODE SYSTEM
+        /// </summary>
+        /// <param name="y0">initial vector</param>
+        /// <param name="start">start time</param>
+        /// <param name="end">end time</param>
+        /// <param name="N">number of subintervals</param>
+        /// <param name="f">ode function</param>
+        /// <returns>approximations</returns>
         public static Vector<double>[] FourthOrder(Vector<double> y0, double start, double end, int N, Func<double, Vector<double>, Vector<double>> 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<double, double, double> ode = (t, y) => t + 2 * y * t;
             Func<double, double> 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
         }
 
         /// <summary>
-        /// Runge-Kutta second order method for first order ODE.
+        /// Runge-Kutta fourth order method for first order ODE.
         /// </summary>
         [Test]
         public void RK4Test()
         {
             Func<double, double, double> ode = (t, y) => t + 2 * y * t;
             Func<double, double> 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