Browse Source

add comments and check the error convergence

pull/394/head
Yoonku Hwang 10 years ago
parent
commit
b2c3fcaf32
  1. 39
      src/Numerics/OdeSolvers/OdeSolvers.cs
  2. 22
      src/UnitTests/OdeSolvers/OdeSolverTest.cs

39
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);

22
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
}
}
}
Loading…
Cancel
Save