I am working with a double integral, and the inner integral has variable bounds. I wrote a function, using SciPy's quad integration, that allows me to evaluate this integral. However, what I want is to just evaluate the inner integral, so that all I am left with is a single, non-evaluated integral with respect to some variable. I want to then plot this "half-way" evaluated double integral vs. a range of that variable so that I can see some trend. However, when I input an array of that variable (it is just 0-10000 but increments of 1, it provides the following error message:
"ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"
I have defined functions before that allow me to input an array (which outputs that function at however many points were in the array), so I am unsure as to why this message is popping up now. I think it has something to do with me using SciPy's "quad" integration when defining the function. How do I go around this so that I can input an array?
import numpy as np
from scipy import integrate
from scipy.integrate import quad
import matplotlib.pyplot as plt
#This is the array of values of the variable I ultimately want to have the function plotted
against
timearray = np.arange(0,10000,1)
#This below defines the general function, that is with respect to two variables (p and t)
def thomtest(p,t):
d = 3.086e22
c = 2.998e10
return (3*((np.cos(p + (2*t/(d*p))))**2))/(8*(t+((p**2)*d/(2*c))))
#The function below evaluates just the inner-integral
def phib(t):
d = 3.086e22
c = 2.998e10
return quad(thomtest,0.00001*(c*t)/d,np.pi, args=(t))[0]
#This evaluates the outer-integral, giving me the complete numerical answer
doubleintegral = quad(phib,0,((3.086e22)/(2.998e10)))
#This below is what gives me the error message: "ValueError: The truth
#value of an array with more than one element is ambiguous.
#Use a.any() or a.all()". Apparently I cannot input an array
print(phib(timearray))
In [2]: doubleintegral
Out[2]: (0.9892936902920587, 1.3643899787751934e-08)
In [3]: phib(0)
Out[3]: 6.377997354641736e-10
In [4]: phib(np.arange(0,5))
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-4-89b7b858c2f5> in <module>
----> 1 phib(np.arange(0,5))
<ipython-input-1-4c64e69e4f62> in phib(t)
10 d = 3.086e22
11 c = 2.998e10
---> 12 return quad(thomtest,0.00001*(c*t)/d,np.pi, args=(t))[0]
13
14 doubleintegral = quad(phib,0,((3.086e22)/(2.998e10)))
/usr/local/lib/python3.6/dist-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
336
337 # check the limits of integration: \int_a^b, expect a < b
--> 338 flip, a, b = b < a, min(a, b), max(a, b)
339
340 if weight is None:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
In [5]:
The error comes when it test the integration bounds, making sure that the first is smaller than the second.
But with an array t
, that test is multivalued, and cannot be used:
In [6]: d = 3.086e22
...: c = 2.998e10
In [7]: 0.00001*(c*np.arange(5))/d
Out[7]:
array([0.00000000e+00, 9.71484122e-18, 1.94296824e-17, 2.91445237e-17,
3.88593649e-17])
With quad
you have to run your code, phib
and doubleintegral
once for each element of t
. In numpy
we try to avoid iterations, but quad
only works with scalar bounds. So good old iteration is required.