Numerical evaluation of line integral

Question

I have written a function that evaluates a contour integral in Python (Poincaré - Hopf index for 2D vector fields). It works, however is really slow if the arg res=N is chosen higher than 50. Input are two arrays of dimension N x N. As the integrals should lead to an integer, (usually -1, 0, or 1) only small accuracy is necessary. I believe that most time is lost doing the interpolation, but I don't know any other method.

I have added a sample vector field \$ \mathbf{M} = (mx, my)^T \$. This produces output -1 if \$ (0, 0) \$ is included in the contour and 0 else. This works as a sample input. Real data is usually less smooth.

import numpy as np
from scipy.interpolate import interp2d
from scipy.integrate import quad
from scipy.misc import derivative

N = 100
xvec = np.linspace(-5, 5, N)
X, Y = np.meshgrid(xvec, xvec)
mx = Y / (X**2 + Y**2)**(1/2)
my = X / (X**2 + Y**2)**(1/2)

def partial_derivative(func, var=0, point=[]):
    args = point[:]

    def wraps(x):
        args[var] = x
        return func(*args)

    return derivative(wraps, point[var], dx=1e-6)


def poin_int(arr_x, arr_y, alpha=5, res=N, r=1, x0=0, y0=0):
    xvec = np.linspace(-alpha, alpha, res)
    X, Y = np.meshgrid(xvec, xvec)  # create grid
    x_func = interp2d(X, Y, arr_x, kind='cubic')  # interpolate data
    y_func = interp2d(X, Y, arr_y, kind='cubic')  # interpolate data

    def _x(t):  # x-coordinate of circular contour
        return np.cos(t) * r + x0

    def _dx(t):  # derivate of _x(t)
        return -np.sin(t) * r

    def _y(t):  # y-coordinate of circular contour
        return np.sin(t) * r + y0

    def _dy(t):  # derivative of _y(t)
        return np.cos(t) * r

    def _integrand1(t):  # integrand of first integral
        return (x_func(_x(t), _y(t)) *
                (partial_derivative(y_func, 0, [_x(t), _y(t)]) * _dx(t) +
                 partial_derivative(y_func, 1, [_x(t), _y(t)]) * _dy(t)) /
                (x_func(_x(t), _y(t))**2 + y_func(_x(t), _y(t))**2)  * 2 * np.pi)

    def _integrand2(t):  # integrand of second integral
        return (-y_func(_x(t), _y(t)) *
                (partial_derivative(x_func, 0, [_x(t), _y(t)]) * _dx(t) +
                 partial_derivative(x_func, 1, [_x(t), _y(t)]) * _dy(t)) /
                (x_func(_x(t), _y(t))**2 + y_func(_x(t), _y(t))**2 * 2 * np.pi))

    _int1 = quad(_integrand1, 0, 2*np.pi, epsabs=1e-2)  
    _int2 = quad(_integrand2, 0, 2*np.pi, epsabs=1e-2)
    return (_int1[0] + _int2[0])  # should be an integer

poin_int(mx, my, res=N, r=1, x0=0, y0=0)

Answer 1

Accuracy

the integrals should lead to an integer, (usually -1, 0, or 1)

Well... they don't. I don't know why this is. Currently the sample input always produces -20.45.

Are you absolutely sure that _integrand2 should have 2*np.pi at a different level of nesting than the same term in _integrand1? That looks like a potential accident.

Also a potential accident is

mx = Y / (X**2 + Y**2)**(1/2)
my = X / (X**2 + Y**2)**(1/2)

Are those really meant to be swapped?

Normally, you should use np.sqrt instead of **(1/2); but since you also have **2 on the inside, use hypot().

Python

Move all of these:

N = 100
xvec = np.linspace(-5, 5, N)
X, Y = np.meshgrid(xvec, xvec)
mx = Y / (X**2 + Y**2)**(1/2)
my = X / (X**2 + Y**2)**(1/2)

out of the global namespace.

Add typehints.

Wrap the upper call to poin_int in a __main__ guard.

The contents of _integrand1 and _integrand2 are far too dense - they need to be split onto multiple lines and multiple intermediate variables.

Deprecations

derivative() is deprecated and about to be removed from Scipy. But also: because your data are already interpolated, it's not exactly the most appropriate method to call. Instead, consider calling into partial_derivative.

interp2d is also deprecated. Normally I'd suggest RegularGridInterpolator, but because you need partial derivatives, use RectBivariateSpline.

Suggested

import numpy as np
from scipy.interpolate import RectBivariateSpline
from scipy.integrate import quad


def poin_int(
    arr_x: np.ndarray,
    arr_y: np.ndarray,
    res: int,
    alpha: float = 5.,
    r: float = 1.,
    x0: float = 0.,
    y0: float = 0.,
) -> float:
    xvec = np.linspace(start=-alpha, stop=alpha, num=res)
    x_func = RectBivariateSpline(xvec, xvec, arr_x)
    y_func = RectBivariateSpline(xvec, xvec, arr_y)
    dx0 = x_func.partial_derivative(dx=1, dy=0)
    dx1 = x_func.partial_derivative(dx=0, dy=1)
    dy0 = y_func.partial_derivative(dx=1, dy=0)
    dy1 = y_func.partial_derivative(dx=0, dy=1)

    def _x(t: float) -> float:  # x-coordinate of circular contour
        return np.cos(t)*r + x0

    def _dx(t: float) -> float:  # derivative of _x(t)
        return -np.sin(t)*r

    def _y(t: float) -> float:  # y-coordinate of circular contour
        return np.sin(t)*r + y0

    def _dy(t: float) -> float:  # derivative of _y(t)
        return np.cos(t)*r

    def integrand1(t: float) -> float:  # integrand of first integral
        p = _x(t), _y(t)
        return (
            x_func(*p)
            * (dy0(*p)*_dx(t) + dy1(*p)*_dy(t))
            / (
                x_func(*p)**2 + y_func(*p)**2
            )
            *2*np.pi
        )

    def integrand2(t: float) -> float:  # integrand of second integral
        p = _x(t), _y(t)
        return (
            -y_func(*p)
            * (dx0(*p)*_dx(t) + dx1(*p)*_dy(t))
            / (
                x_func(*p)**2 + y_func(*p)**2 * 2*np.pi  # really?
            )
        )

    int1, err = quad(func=integrand1, a=0, b=2*np.pi, epsabs=1e-4)
    int2, err = quad(func=integrand2, a=0, b=2*np.pi, epsabs=1e-4)
    return int1 + int2  # should be an integer


def demo() -> None:
    N = 100
    xvec = np.linspace(start=-5, stop=5, num=N)
    X, Y = np.meshgrid(xvec, xvec)
    denom = np.hypot(X, Y)
    mx = X/denom
    my = Y/denom
    result = poin_int(arr_x=mx, arr_y=my, res=N, r=1, x0=0, y0=0)
    print(result)


if __name__ == '__main__':
    demo()