sebastiaagramunt / micromagnetics Goto Github PK

def derivative_coefficients(stencil: List[float], deg: int):
    """
    Calculates the coefficients for the numerical derivative
    Inputs:
        deg: int -> the max degree of the derivatives
    Output:
        The coefficients to calculate the derivative at one point. 
    Refs:
        http://web.media.mit.edu/~crtaylor/calculator.html
        https://en.wikipedia.org/wiki/Finite_difference_coefficient
    """

    # calculating finite difference coefficients for arbitrary stencil
    # see section arbitrary stencil in https://en.wikipedia.org/wiki/Finite_difference_coefficient
    N = len(stencil)
    stencil_matrix = np.zeros((N, N), dtype=np.int64)
    for i, j in [(x, y) for x in range(N) for y in range(N)]:
        stencil_matrix[i][j] = pow(stencil[j], i)

    stencil_matrix_inv = np.linalg.inv(stencil_matrix)
    v = np.array(
        [0 if i != deg else np.math.factorial(deg) for i in range(N)], dtype=np.int
    )

    return stencil_matrix_inv.dot(v)

def numeric_derivative(target_fn: Callable, stencil: List[float], a: float, h: float, deg: int):

    """
    Calculate numerical d derivative at point a using points in self.strencil
    Input:
        a: point where we want to calculate the derivative
        h: separation between evaluation points
        deg: degree of the numeric derivative
    """

    coeffs = derivative_coefficients(deg)
    res = 0
    for s, coef in zip(stencil, coeffs):
        res += coef * target_fn(a + s * h) / pow(h, deg)

    return res