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
    It would be very helpful to include some input to the functions. That way, debugging and optimization is possible. – maxb Aug 14 at 5:56

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.