Browse Source

done Adams-Bashforth second order method to solve ODE

pull/397/head
Yoonku Hwang 10 years ago
parent
commit
7de883a41a
  1. 14
      src/Numerics/OdeSolvers/AdamsBashforth.cs
  2. 2
      src/UnitTests/OdeSolvers/OdeSolverTest.cs

14
src/Numerics/OdeSolvers/AdamsBashforth.cs

@ -55,6 +55,7 @@ namespace MathNet.Numerics.OdeSolvers
}
return y;
}
/// <summary>
/// Second Order AB Method(Require two initial guesses)
/// </summary>
@ -64,23 +65,22 @@ namespace MathNet.Numerics.OdeSolvers
/// <param name="N">Number of subintervals</param>
/// <param name="f">ode model</param>
/// <returns></returns>
public static double[] SecondOrder(double y0, double start, double end, int N, Func<double, double, double> f)
public static double[] SecondOrder(double y0, double start, double end, int N, Func<double, double, double> f)
{
double dt = (end - start) / (N - 1);
double t = start;
double[] y = new double[N];
double fold = f(t, y0);
double ytilde = y0 + (2*dt/3) * fold;
double y1 = y0 + (dt/4) * fold + (3*dt/4) * f(t + (2*dt/3), ytilde );
double fold = f(t, y0);
double ytilde = y0 + 2.0 * dt / 3.0 * fold;
double y1 = y0 + dt / 4.0 * fold + 3.0 * dt / 4 * f(t + 2.0 / 3.0 * dt, ytilde);
y[0] = y0;
t += dt;
y[1] = y1;
for (int i = 2; i < N; i++)
{
y[i] = y1 + dt * (1.5 * f(t + dt, y1) - 0.5 * f(t, y0));
t += dt;
y0 = y[i-1];
y0 = y[i - 1];
y1 = y[i];
}
return y;

2
src/UnitTests/OdeSolvers/OdeSolverTest.cs

@ -72,7 +72,7 @@ namespace MathNet.Numerics.UnitTests.OdeSolvers
for (int k = 0; k < 4; k++)
{
double y0 = 0;
double[] y_t = AdamsBashforth.SecondOrder(y0, 0, 2, Convert.ToInt32(Math.Pow(2, k + 6)), ode);
double[] y_t = AdamsBashforth.SecondOrder(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);

Loading…
Cancel
Save