From 02a0af61a236ebf68a47e6b15d0d33945a7a2bd7 Mon Sep 17 00:00:00 2001
From: Yoonku Hwang <yz0624@gmail.com>
Date: Fri, 13 May 2016 20:23:46 +0900
Subject: [PATCH] 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<double, double, double> 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()));
             }
         }
+
+        /// <summary>
+        /// Runge-Kutta second 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);
+            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