# Numerical evaluation of line integral

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.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)

• It would be very helpful to include some input to the functions. That way, debugging and optimization is possible. – maxb Aug 14 '18 at 5:56