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