quadpy
quadpy copied to clipboard
sigma-py
Double integration issue - passing numpy array instead of single values
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.
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.
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
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.
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]
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.
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
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):?
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.
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.
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?
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.
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.