Did I properly implement double integral equation in dblquad Python? Helicoidal Surface Theory [closed]

I show the equations (simplified Helicoidal Surface Theory for implementation purposes) that I want to calculate numerically using Python.

and the code:

from scipy.integrate import quad
import numpy as np
from scipy import interpolate
import time
start_time = time.time()

input="-0.5 0.0 \
-0.3 0.9    \
0.0 0.8 \
0.3 0.4 \
0.5 0.02"

input_coordinates = np.genfromtxt(input.splitlines()).reshape(-1,2) # shape to 2 columns, any number of rows
x_coordinates = input_coordinates[:,0]
H_values = input_coordinates[:,1]
H_interpolation = interpolate.InterpolatedUnivariateSpline(x_coordinates, H_values)

def complex_dblquad(func, a, b, g, h, **kwargs):
def real_func(z, x):
return np.real(func(z, x))
def imag_func(z, x):
return np.imag(func(z, x))
real_integral = dblquad(real_func, a, b, g, h, **kwargs)
imag_integral = dblquad(imag_func, a, b, g, h, **kwargs)
return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:])

complex_integral = complex_dblquad(lambda z,x: np.sqrt(1+z*z)**2*(2/np.sqrt(1+z*z))**2*H_interpolation(x)*np.exp(1j*2/np.sqrt(1+z*z)*x), 0, 1, -0.5, 0.5)

print("--- %s seconds ---" % (time.time() - start_time))


The question is - is it implemented correctly? If not - what is wrong?

• I don't think you should be dumping your code here to get it checked for correctness. If you think there is a mistake then say why you think there is a mistake, or focus on one specific, short section of code. Jun 9, 2021 at 14:09

It looks correct, but it is difficult to interpret it is not a good practice to expand all the variables in one single expression, you can use the use definitions that match your equations more easily.

def integrand_function(z, x):
Mr = np.sqrt(1 + z**2)
kx = 2/Mr;
# merging the integrand of Psi in order to use the dblquad function.
return 2 * Mr**2 * kx**2 H_interpolation(x)*np.exp(1j * kx * x)


Then your complex integral is evaluated with

complex_integral = complex_dblquad(integrand_function, 0, 1, -0.5, 0.5)


And hopefully you can prove more easily the correctness this way.

• Thank you very much! It helps a lot! Jun 8, 2021 at 11:34