everyone. I was trying to compute an integral numerically using scipy.quad. However, it returns a really huge array with no other error. Can anyone see the problem? The code works very fine with simple integral like quad(lambda x:x**2,0,4) but for the code below try a=20,b=20,c=20,j=1 you get a really huge list of result with repeating numbers
def Lj(a,b,c,j):
R1=a
R2=b
R3=c
if j==1:
Rj=a
elif j==2:
Rj=b
elif j==3:
Rj=c
else:
raise ValueError('j needs to be an integer between 1 and 3')
result = (R1*R2*R3)*quad(lambda s: 1/((s+Rj**2)*sqrt(((s+R1**2)*(s+R2**2)*(s+R3**2))))/2, 0, inf)
return result
quad
returns a tuple containing the integral and an estimate of the error. You multiply this tuple by R1*R2*R3
, which replicates the tuple R1*R2*R3
times.
To fix this, use just the first value returned by quad
. That is, replace this:
result = (R1*R2*R3)*quad(lambda s: 1/((s+Rj**2)*sqrt(((s+R1**2)*(s+R2**2)*(s+R3**2))))/2, 0, inf)
with something like this:
intgrl, err = quad(lambda s: 1/((s+Rj**2)*sqrt(((s+R1**2)*(s+R2**2)*(s+R3**2))))/2, 0, inf)
result = (R1*R2*R3)*intgrl