bdf2 adaptive

def emr_step(dt_n, y_n, dy_n, dt_nm1, y_nm1):
    """Take a single step of the explicit midpoint rule.
    From G&S pg. 715 and Prinja's thesis pg.45.
    """
    dtr = dt_n / dt_nm1
    y_np1 = (1 - dtr**2)*y_n + (1 + dtr)*dt_n*dy_n + (dtr**2)*(y_nm1)
    return y_np1

def bdf2_dydt(ts, ys):
    """Get dy/dt at time ts[-1] (allowing for varying dt).
    Gresho & Sani, pg. 715"""
    dt_n = ts[-1] - ts[-2]
    dt_nm1 = ts[-2] - ts[-3]

    y_np1 = ys[-1]
    y_n = ys[-2]
    y_nm1 = ys[-3]

    # Copied from oomph-lib (algebraic rearrangement of G&S forumla).
    weights = [1.0/dt_n + 1.0/(dt_n + dt_nm1),
               - (dt_n + dt_nm1)/(dt_n * dt_nm1),
               dt_n / ((dt_n + dt_nm1) * dt_nm1)]

    dydt = sp.dot(weights, [y_np1, y_n, y_nm1])

    return dydt


def bdf2_mp_gs_lte_estimate(ts, ys):
    """Estimate LTE using combination of bdf2 and explicit midpoint. From
    oomph-lib and G&S.
    """

    # Get local values (makes maths more readable)
    dt_n = ts[-1] - ts[-2]
    dt_nm1 = ts[-2] - ts[-3]
    dtr = dt_n / dt_nm1
    dtrinv = 1.0 / dtr

    y_np1 = ys[-1]
    y_n = ys[-2]
    y_nm1 = ys[-3]

    # Invert bdf2 to get predictor data (using the exact same function as
    # was used in the residual calculation).
    dy_n = bdf2_dydt(ts, ys)

    # Calculate predictor value (variable dt explicit mid point rule)
    y_np1_EMR = emr_step(dt_n, y_n, dy_n, dt_nm1, y_nm1)

    error_weight = ((1.0 + dtrinv)**2) / \
        (1.0 + 3.0*dtrinv + 4.0 * dtrinv**2
         + 2.0 * dtrinv**3)

    # Calculate truncation error -- oomph-lib
    error = (y_np1 - y_np1_EMR) * error_weight

    return error

emr_step corresponds to the calculations done by BDF<2>::calculate_predicted_values with the weights set by BDF<2>::set_predictor_weights().

bdf2_dydt corresponds to TimeStepper::time_derivative(1,...) with the weights set by BDF<2>::set_weights().

bdf2_mp_gs_lte_estimate corresponds to BDF<2>::temporal_error_in_value with the error weight set by BDF<2>::set_error_weights().

In python -1 gets the last element of the list, -2 the one before last, etc.

I store time and y values in time order, i.e. ts[0] = t_0, ts[-2] = t_n and ts[-1] = t_{n+1}.

** is the power operator