Math.py 文件源码

def fit_quartic(y0, y1, g0, g1):
    """Fit constrained quartic polynomial to function values and erivatives at x = 0,1.

    Returns position and function value of minimum or None if fit fails or has
    a maximum. Quartic polynomial is constrained such that it's 2nd derivative
    is zero at just one point. This ensures that it has just one local
    extremum.  No such or two such quartic polynomials always exist. From the
    two, the one with lower minimum is chosen.
    """
    def g(y0, y1, g0, g1, c):
        a = c+3*(y0-y1)+2*g0+g1
        b = -2*c-4*(y0-y1)-3*g0-g1
        return np.array([a, b, c, g0, y0])

    def quart_min(p):
        r = np.roots(np.polyder(p))
        is_real = np.isreal(r)
        if is_real.sum() == 1:
            minim = r[is_real][0].real
        else:
            minim = r[(r == max(-abs(r))) | r == -max(-abs(r))][0].real
        return minim, np.polyval(p, minim)

    D = -(g0+g1)**2-2*g0*g1+6*(y1-y0)*(g0+g1)-6*(y1-y0)**2  # discriminant of d^2y/dx^2=0
    if D < 1e-11:
        return None, None
    else:
        m = -5*g0-g1-6*y0+6*y1
        p1 = g(y0, y1, g0, g1, .5*(m+np.sqrt(2*D)))
        p2 = g(y0, y1, g0, g1, .5*(m-np.sqrt(2*D)))
        if p1[0] < 0 and p2[0] < 0:
            return None, None
        [minim1, minval1] = quart_min(p1)
        [minim2, minval2] = quart_min(p2)
        if minval1 < minval2:
            return minim1, minval1
        else:
            return minim2, minval2