Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Double integration issue - passing numpy array instead of single values #414

Open
tom634 opened this issue May 20, 2021 · 16 comments
Open

Comments

@tom634
Copy link

tom634 commented May 20, 2021

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.

@tom634 tom634 changed the title Double integration Double integration issue - passing numpy array instead of single values May 20, 2021
@nschloe
Copy link
Collaborator

nschloe commented May 20, 2021

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
Copy link
Author

tom634 commented May 20, 2021

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
Copy link
Collaborator

nschloe commented May 20, 2021

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
Copy link
Author

tom634 commented May 20, 2021

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
Copy link
Collaborator

nschloe commented May 21, 2021

It's scipy.quad that gives that error.

@tom634
Copy link
Author

tom634 commented May 21, 2021

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
Copy link
Collaborator

nschloe commented May 21, 2021

You can loop then.

@tom634
Copy link
Author

tom634 commented May 21, 2021

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
Copy link
Collaborator

nschloe commented May 21, 2021

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

@tom634
Copy link
Author

tom634 commented May 21, 2021

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
Copy link
Collaborator

nschloe commented May 21, 2021

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
Copy link
Author

tom634 commented May 21, 2021

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
Copy link
Author

tom634 commented Jun 2, 2021

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
Copy link
Author

tom634 commented Jun 7, 2021

In order to simplify the above question, I show the equations (simplified for implementation purposes) that I want to calculate numerically.
obraz
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.

@tom634
Copy link
Author

tom634 commented Jul 22, 2021

Hi @nschloe do you know maybe how to resolve this issue?

@nschloe
Copy link
Collaborator

nschloe commented Jul 22, 2021

I don't have much time these days. I can point you to https://github.com/nschloe/quadpy/wiki/Dimensionality-of-input-and-output-arrays.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants