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.

Update: 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)
  • 1
    \$\begingroup\$ It would be very helpful to include some input to the functions. That way, debugging and optimization is possible. \$\endgroup\$ – maxb Aug 14 '18 at 5:56

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.