Browse Source

blowing up error when adding second order AB method

pull/397/head
Yoonku Hwang 10 years ago
parent
commit
75dd96c987
  1. 34
      src/Numerics/OdeSolvers/AdamsBashforth.cs
  2. 26
      src/UnitTests/OdeSolvers/OdeSolverTest.cs

34
src/Numerics/OdeSolvers/AdamsBashforth.cs

@ -27,15 +27,11 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace MathNet.Numerics.OdeSolvers
{
public static class AdamsBashforth
{
/// <summary>
/// First Order AB method(same as Forward Euler)
/// </summary>
@ -59,10 +55,32 @@ namespace MathNet.Numerics.OdeSolvers
}
return y;
}
public static double[] SecondOrder()
/// <summary>
/// Second Order AB Method(Require two initial guesses)
/// </summary>
/// <param name="y0">Initial value 1</param>
/// <param name="y1">Initial value 2</param>
/// <param name="start">Start Time</param>
/// <param name="end">End Time</param>
/// <param name="N">Number of subintervals</param>
/// <param name="f">ode model</param>
/// <returns></returns>
public static double[] SecondOrder(double y0, double y1, 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];
y[0] = y0;
y[1] = y1;
for (int i = 2; i < N; i++)
{
y1 = f(t + dt, y1);
y[i] = y1 + dt * (1.5 * f(t + dt, y1) - 0.5 * f(t, y0));
t += dt;
y0 = y1;
y1 = y[i];
}
return y;
}
public static double[] ThirdOrder()
@ -75,4 +93,4 @@ namespace MathNet.Numerics.OdeSolvers
return null;
}
}
}
}

26
src/UnitTests/OdeSolvers/OdeSolverTest.cs

@ -27,9 +27,9 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using MathNet.Numerics.OdeSolvers;
using NUnit.Framework;
using System;
using MathNet.Numerics.OdeSolvers;
using System.Linq;
namespace MathNet.Numerics.UnitTests.OdeSolvers
@ -60,6 +60,30 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers
}
Assert.AreEqual(1, ratio, 0.01);// Check error convergence order
}
[Test]
public void AB2Test()
{
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;
int N = Convert.ToInt32(Math.Pow(2, k + 6));
double dt = 2.0 / (N - 1);
double[] y_t = AdamsBashforth.SecondOrder(y0, sol(dt), 0, 2, N, 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(2, ratio, 0.01);// Check error convergence order
}
/// <summary>
/// Runge-Kutta second order method for first order ODE.
/// </summary>

Loading…
Cancel
Save