Double integration issue - passing numpy array instead of single values

[tom634]

I want to calculate integral located internally of the integral with imaginary numbers.

My code:

from scipy.integrate import quad
import numpy as np
import scipy
from scipy.special import jv
import math
import quadpy

def M_r(z):
    print("z:",z)
    return math.sqrt(z*z)

def integral_P(z):
    print("Calling function, z=",z)
    return (M_r(z)*np.exp(1j))

Integral = quadpy.quad(integral_P, 0, 1, limit=10)

print("Quadpy",Integral)

When I use scipy quad integration:
Integral = quad(integral_P, 0, 1, limit=10)
In the printout z values are single values, nor numpy array:

Calling function, z= 0.013046735741414128
z: 0.013046735741414128
Calling function, z= 0.9869532642585859
z: 0.9869532642585859
Calling function, z= 0.06746831665550773
z: 0.06746831665550773

but on other hand I have casting problem:

ComplexWarning: Casting complex values to real discards the imaginary part
  return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)

So I want to use quadpy.quad to deal with imaginary numbers. When I run this code, I have error:

Exception has occurred: TypeError
only size-1 arrays can be converted to Python scalars

And the printout is:

Calling function, z= [0.00217142 0.01304674 0.03492125 0.06746832 0.10959114 0.16029522
 0.21862143 0.2833023  0.35280357 0.42556283 0.5        0.57443717
 0.64719643 0.7166977  0.78137857 0.83970478 0.89040886 0.93253168
 0.96507875 0.98695326 0.99782858]
z: [0.00217142 0.01304674 0.03492125 0.06746832 0.10959114 0.16029522
 0.21862143 0.2833023  0.35280357 0.42556283 0.5        0.57443717
 0.64719643 0.7166977  0.78137857 0.83970478 0.89040886 0.93253168
 0.96507875 0.98695326 0.99782858]

Quadpy should pass single z values, over iterations, not entire numpy array. How can I solve this issue? The type of z is:
z: <class 'numpy.ndarray'>

After some tests I noticed that numpy array is passed also to internal M_r(z) function even if in the main function integral_P there are no imaginary numbers.

[nschloe]

For me, your code correctly gives

[0.5  0.25]
Calling function, z= [0.00217142 0.01304674 0.03492125 0.06746832 0.10959114 0.16029522
 0.21862143 0.2833023  0.35280357 0.42556283 0.5        0.57443717
 0.64719643 0.7166977  0.78137857 0.83970478 0.89040886 0.93253168
 0.96507875 0.98695326 0.99782858]
z: [0.00217142 0.01304674 0.03492125 0.06746832 0.10959114 0.16029522
 0.21862143 0.2833023  0.35280357 0.42556283 0.5        0.57443717
 0.64719643 0.7166977  0.78137857 0.83970478 0.89040886 0.93253168
 0.96507875 0.98695326 0.99782858]

and then errors out because math.sqrt cannot handle vectors.

[tom634]

I changed math.sqrt to np.sqrt, but now I have another problem - quad divided to real and imaginary parts and quadpy calculating both parts give different results - why?
Here is code:

from scipy.integrate import quad
import numpy as np
import scipy
from scipy.special import jv
import math
import quadpy


B=8
m=4
M_x = 0.3
M_T = 1
theta=0.3

def M_r(z):
    print("M_r",z)
    return np.sqrt(z*z)

def phi_0(z):
    print("phi:",z)
    return 2*m*B/(z*M_r(z))

def integral_P_Vm(z):
    print("Calling function, z=",z)
    # return 2*np.exp(1j*(z))*jv(m*B, (B*z))*2
    return M_r(z)**2*np.exp(1j*(phi_0(z)))*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2
    # return np.sin(z) # test function

def integral_P_Vm_real(z):
    print("Calling real function, z=",z)
    # return np.real(2*np.exp(1j*(z))*jv(m*B, (B*z))*2)
    return np.real(M_r(z)**2*np.exp(1j*(phi_0(z)))*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2)

def integral_P_Vm_imag(z):
    print("Calling imag function, z=",z)
    # return np.imag(2*np.exp(1j*(z))*jv(m*B, (B*z))*2)
    return np.imag(M_r(z)**2*np.exp(1j*(phi_0(z)))*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2)

P_Vm_integral = quadpy.quad(integral_P_Vm, 0, 1, limit=100)
P_Vm_integral_real = quad(integral_P_Vm_real, 0, 1, limit=100)
P_Vm_integral_imag = quad(integral_P_Vm_imag, 0, 1, limit=100)


print("Quadpy",P_Vm_integral)
print("Quad:", P_Vm_integral_real[0]+1j*P_Vm_integral_imag[0], P_Vm_integral_real[1]+P_Vm_integral_imag[1])

The output is:

Quadpy ((-8.875555588362743e-12-7.735411867050377e-13j), 1.215938850446882e-11)
Quad: (-1.833732372097316e-12+1.7547331112662863e-12j) 1.8061540068960578e-11

I observed that if I change phi_0 to:

def phi_0(z):
    print("phi:",z)
    return 2/(z*M_r(z))

results are almost same:

Quadpy ((-6.092835777513197e-12+9.738512507900862e-12j), 9.955172432146452e-17)
Quad: (-6.0928069104517275e-12+9.738517119655075e-12j) 1.4288475456681177e-16

[nschloe]

Different algorithms, different results. Both results are in the order of machine precision, and they differ by less than 10^{-11}. The error indicates that, too. Everything seems in order.

[tom634]

Thank you for the comment. I have another issue with integrating integral. My code is:

from scipy.integrate import quad
import numpy as np
import scipy
from scipy.special import jv
import math
import quadpy

B=8
m=4
M_x = 0.3
M_T = 1
theta=0.3

def f_D(x, z):
    print("f_D:",z)
    return 1*np.exp(1j*x)*z

def M_r(z):
    print("M_r",z)
    return np.sqrt(z*z)

def psi_D(z):
    print("psi_D",z)
    return quad(f_D, -0.5, 0.5, limit=50, args=(z))[0]

def integral_P_Vm(z):
    print("Calling function, z=",z)
    return M_r(z)**2*np.exp(1j)*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2*psi_D(z)

Quad_integral = quad(integral_P_Vm, 0, 1, limit=100)
Quadpy_integral = quadpy.quad(integral_P_Vm, 0, 1, limit=100)

print("Quad",Quad_integral)
print("Quadpy",Quadpy_integral)

Scipy Quad integrates (but with no imaginary parts).
Integrating with Quadpy gives error in line number 22:
return quad(f_D, -0.5, 0.5, limit=50, args=(z))[0]
Error:

Exception has occurred: TypeError
only size-1 arrays can be converted to Python scalars

The last output is following:

psi_D [0.00217142 0.01304674 0.03492125 0.06746832 0.10959114 0.16029522
 0.21862143 0.2833023  0.35280357 0.42556283 0.5        0.57443717
 0.64719643 0.7166977  0.78137857 0.83970478 0.89040886 0.93253168
 0.96507875 0.98695326 0.99782858]
f_D: [0.00217142 0.01304674 0.03492125 0.06746832 0.10959114 0.16029522
 0.21862143 0.2833023  0.35280357 0.42556283 0.5        0.57443717
 0.64719643 0.7166977  0.78137857 0.83970478 0.89040886 0.93253168
 0.96507875 0.98695326 0.99782858]

[nschloe]

It's scipy.quad that gives that error.

[tom634]

When I change:
return quad(f_D, -0.5, 0.5, limit=50, args=(z))[0]
to
return quadpy.quad(f_D, -0.5, 0.5, limit=50, args=(z))[0]
I have error:

Exception has occurred: TypeError
f_D() takes 2 positional arguments but 22 were given

the
print("psi_D",z)
command gives the output:

psi_D [0.00217142 0.01304674 0.03492125 0.06746832 0.10959114 0.16029522
 0.21862143 0.2833023  0.35280357 0.42556283 0.5        0.57443717
 0.64719643 0.7166977  0.78137857 0.83970478 0.89040886 0.93253168
 0.96507875 0.98695326 0.99782858]

The issue is that single values of z should be passed to def f_D(x, z): during integrals.

[nschloe]

You can loop then.

[tom634]

I tried looping with an error:
Exception has occurred: TypeError f_D() argument after * must be an iterable, not numpy.float64

def psi_D(z):
    print("psi_D",z)
    for e in z:
        print(e)
        return quadpy.quad(f_D, -0.5, 0.5, limit=50, args=(e))[0]

Output is:

psi_D [0.00217142 0.01304674 0.03492125 0.06746832 0.10959114 0.16029522
 0.21862143 0.2833023  0.35280357 0.42556283 0.5        0.57443717
 0.64719643 0.7166977  0.78137857 0.83970478 0.89040886 0.93253168
 0.96507875 0.98695326 0.99782858]
0.0021714184870961772

[nschloe]

Always post the full code, I don't want to make any guesses.

[tom634]

Full code:

from scipy.integrate import quad
import numpy as np
import scipy
from scipy.special import jv
import math
import quadpy

B=8
m=4
M_x = 0.3
M_T = 1
theta=0.3

def f_D(x, *args):
    print("f_D:",args)
    return 1*np.exp(1j*x)*z

def M_r(z):
    print("M_r",z)
    return np.sqrt(z*z)

def psi_D(z):
    print("psi_D",z)
    for e in z:
        print(e)
        return quadpy.quad(f_D, -0.5, 0.5, limit=50, args=(e))[0]

def integral_P_Vm(z):
    print("Calling function, z=",z)
    return M_r(z)**2*np.exp(1j)*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2*psi_D(z)

# Quad_integral = quad(integral_P_Vm, 0, 1, limit=100)
Quadpy_integral = quadpy.quad(integral_P_Vm, 0, 1, limit=100)

# print("Quad",Quad_integral)
print("Quadpy",Quadpy_integral)

Error is at line 26:
return quadpy.quad(f_D, -0.5, 0.5, limit=50, args=(e))[0]
Error code is:

Exception has occurred: TypeError
f_D() argument after * must be an iterable, not numpy.float64

Also I am not sure if the integral
return quadpy.quad(f_D, -0.5, 0.5, limit=50, args=(e))[0]
is performed only for one z value and then will go back to
def integral_P_Vm(z):?

[nschloe]

What is

    for e in z:
        print(e)
        return quadpy.quad(f_D, -0.5, 0.5, limit=50, args=(e))[0]

? You return after the first iteration.

FYI: I won't have any more time to look at this problem until next week or so.

[tom634]

This is loop over z values in order to calculate internal function f_D(x, z) in order to calculate integral of def psi_D(z): . But currently I have no idea how to prevent stopping iterations (by returning) after first iteration of z and calculate entire z numpy array.

scipy.quad handles loops automatically, but it doesn't handle imaginary numbers. How should I loop in quadpy manually? Currently I have

Exception has occurred: TypeError
f_D() argument after * must be an iterable, not numpy.float64

error.

FYI: OK.

[tom634]

I solved above problem by changing (z) to [z]:

from scipy.integrate import quad
import numpy as np
import scipy
from scipy.special import jv
import math
import quadpy

B=8
m=4
M_x = 0.3
M_T = 1
theta=0.3

def f_D(x, z):
    print("f_D:")
    return 1*np.exp(1j*x)*z

def M_r(z):
    print("M_r",z)
    return np.sqrt(z*z)

def psi_D(z):
    print("psi_D",z)
    return quadpy.quad(f_D, -0.5, 0.5, limit=50, args=[z])[0]

def integral_P_Vm(z):
    print("Calling function, z=",z)
    return M_r(z)**2*np.exp(1j)*jv(m*B, ((m*B*z*np.sin(theta)/(1-M_x*np.cos(theta)))))*2*psi_D(z)

# Quad_integral = quad(integral_P_Vm, 0, 1, limit=100)
Quadpy_integral = quadpy.quad(integral_P_Vm, 0, 1, limit=100)

# print("Quad",Quad_integral)
print("Quadpy",Quadpy_integral)

But now I have another issue.

I have following code:

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

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 k_x(z, M_r):
    print("k_x(z, M_r)")
    return (2)/(M_r(z))

def H(x, z):
    print("H(x, z)")
    print("X shape:", x.shape, "Current X", x)
    print("Z shape", z.shape, "Current z:", z)
    return H_interpolation(x)*np.exp(1j*k_x(z, M_r))#*x)

def f_D(x, z):
    print("f_D")
    return 1*np.exp(1j*x)*z

def M_r(z):
    print("M_r")
    return np.sqrt(z*z)

def psi_D(z):
    print("psi_D")
    return quadpy.quad(H, -0.5, 0.5, limit=50, args=[z])[0]

def integral_P_Vm(z):
    print("Calling function")
    return M_r(z)**2*np.exp(1j*psi_D(z))

Quadpy_integral = quadpy.quad(integral_P_Vm, 0, 1, limit=100)

print("Quadpy",Quadpy_integral)

At first I have shape (21,) for both variables x and z. But during calculations the shape of x variable grows to (42,) and I have error:

Exception has occurred: ValueError
operands could not be broadcast together with shapes (42,) (21,) 

in the line:
return H_interpolation(x)*np.exp(1j*k_x(z, M_r))#*x)
Why do number of parameters in x variable grow from 21 to 42?

[tom634]

In order to simplify the above question, I show the equations (simplified for implementation purposes) that I want to calculate numerically.

and the code:

from scipy.integrate import quad
import numpy as np
from scipy import interpolate
import quadpy
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 H(x, k_x):
    return H_interpolation(x)*np.exp(1j*k_x)

def psi_V(k_x):
    return quadpy.quad(H, -0.5, 0.5, limit=100, args=[k_x])[0]

def integral_P_Vm(z):
    M_r=np.sqrt(z*z)
    k_x=2/M_r
    return M_r**2*k_x**2*psi_V(k_x)

P_Vm_integral = quadpy.quad(integral_P_Vm, 0, 1, limit=100)

print("Quad",P_Vm_integral)

the above code gives error:
Exception has occurred: ValueError
operands could not be broadcast together with shapes (42,) (21,)
that is because integral
return quadpy.quad(H, -0.5, 0.5, limit=100, args=[k_x])[0]
sometimes is over 21 and sometimes over 42 parameters. How to stabilize number of limit? Is the limit the issue? Changing limit doesn't help.