coefficients: small offsets could raise `Cannot compute accuracy' error unexpectedly

[ajz34]

Hi devs! This repo looks cool :sunglasses:

When I try to calculate finite difference coefficients, I encountered exception `Cannot compute accuracy'. https://github.com/maroba/findiff/blob/e8ca33707e3e25d76bf0f93b2391e466209287b1/findiff/coefs.py#L225

My offsets are somehow small:

import numpy as np
import findiff

offsets = np.array([-4, -2, -1, 0, 1, 2, 4]) * 0.0001
findiff.coefficients(deriv=3, offsets=list(offsets))

The following code may probably alleviate the problem.

def _calc_accuracy(offsets, coefs, deriv, symbolic=False):

    n = deriv + 1
    max_n = 999
    if symbolic:
        break_cond = lambda b: b != 0
        while True:
            b = 0
            for o, coef in zip(offsets, coefs):
                b += coef * o ** n
            if break_cond(b):
                break
            n += 1
            if n > max_n:
                raise Exception('Cannot compute accuracy.')
    else:
        break_cond = lambda b: abs(b) > 1.E-6
        coefs_array = np.array(coefs)
        offsets_scaled = np.array(offsets)
        offsets_scaled = offsets_scaled / np.abs(offsets_scaled).max()
        while True:
            l = coefs_array * offsets_scaled**n
            b = l.sum() / np.abs(l).sum()
            if break_cond(b):
                break
            n += 1
            if n > max_n:
                raise Exception('Cannot compute accuracy.')

    return round(n - deriv)

[maroba]

The offsets argument in the findiff.coefficients function is supposed to be a list of integers, because the offsets are integer multiples of the grid spacing. What are you trying to calculate? Maybe I can point you to the solution how to do that with findiff...

[ajz34]

I'm just tring findiff for fun, not for production solution; and a workaround already exists. So not a big problem to me :joy:

I require offsets like [-1.4h, -h, 0, h, 1.4h] to generate a third derivative. Float 1.4 or derivative order could be any other value. And then I found findiff.coefficients could handle float offsets as large as 1.4 properly, so I tried to feed h=0.0001 into the function, thus got this exception message.

The reason I use float offsets is to generate generalized Rutishauser–Romberg table. Ordinary Rutishauser–Romberg table (or Romberg's Method Wikipedia) only requires offsets like [..., -8h, -4h, -2h, -h, 0, h, 2h, 4h, 8h, ...], and that could be done by findiff. Offsets in the so-called generalized Rutishauser–Romberg table (JMS 2007) is something like [..., -1.4^2 h, -1.4 h, -h, 0, h, 1.4 h, 1.4^2 h, ...]. So I have to handle non-integer offsets.

I only need coefficients, and not considering more advanced features in findiff currently; since the function is not simple as well as really hard to interface with python. So I tried to solve this by myself:

def calculate_findiff_coefs(offsets, deriv):
    offsets = np.asarray(offsets)
    matrix = np.array([offsets**n for n in range(len(offsets))])
    rhs = np.zeros(len(offsets))
    rhs[deriv] = np.math.factorial(deriv)
    return np.linalg.solve(matrix, rhs)

This situation is somehow specialized. For your consideration, I conjucture that probably most audience of this package is to solve differential equations or classroom illustration; non-integer grid is rarely required. I may close this issue.