Browse Source

add Adams-Bashforth 3rd order and verify the order

pull/397/head
Yoonku Hwang 10 years ago
parent
commit
0476e390db
  1. 41
      src/Numerics/OdeSolvers/AdamsBashforth.cs
  2. 21
      src/UnitTests/OdeSolvers/OdeSolverTest.cs

41
src/Numerics/OdeSolvers/AdamsBashforth.cs

@ -57,7 +57,7 @@ namespace MathNet.Numerics.OdeSolvers
}
/// <summary>
/// Second Order AB Method(Require two initial guesses)
/// Second Order AB Method
/// </summary>
/// <param name="y0">Initial value 1</param>
/// <param name="start">Start Time</param>
@ -87,13 +87,46 @@ namespace MathNet.Numerics.OdeSolvers
return y;
}
public static double[] ThirdOrder()
/// <summary>
/// Third Order AB Method
/// </summary>
/// <param name="y0">Initial value 1</param>
/// <param name="start">Start Time</param>
/// <param name="end">End Time</param>
/// <param name="N">Size of output array(the larger, the finer)</param>
/// <param name="f">ode model</param>
/// <returns>approximation with size N</returns>
public static double[] ThirdOrder(double y0, double start, double end, int N, Func<double, double, double> f)
{
return null;
double dt = (end - start) / (N - 1);
double t = start;
double[] y = new double[N];
double k1 = 0;
double k2 = 0;
double k3 = 0;
double k4 = 0;
y[0] = y0;
for (int i = 1; i < 3; i++)
{
k1 = dt * f(t, y0);
k2 = dt * f(t + dt / 2, y0 + k1 / 2);
k3 = dt * f(t + dt / 2, y0 + k2 / 2);
k4 = dt * f(t + dt, y0 + k3);
y[i] = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
t += dt;
y0 = y[i];
}
for (int i = 3; i < N; i++)
{
y[i] = y[i - 1] + dt * (23 * f(t, y[i - 1]) - 16 * f(t - dt, y[i - 2]) + 5 * f(t - 2 * dt, y[i - 3])) / 12.0;
t += dt;
}
return y;
}
/// <summary>
/// Fourth Order AB Method(Require two initial guesses)
/// Fourth Order AB Method
/// </summary>
/// <param name="y0">Initial value 1</param>
/// <param name="start">Start Time</param>

21
src/UnitTests/OdeSolvers/OdeSolverTest.cs

@ -82,6 +82,27 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers
Assert.AreEqual(2, ratio, 0.01);// Check error convergence order
}
[Test]
public void AB3Test()
{
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 = AdamsBashforth.ThirdOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 8)), ode);
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(3, ratio, 0.01);// Check error convergence order
}
[Test]
public void AB4Test()
{

Loading…
Cancel
Save