Reputation: 944
from __future__ import division
import functools
import warnings
import numpy as np
import scipy as sp
from scipy import integrate
from numpy import exp, pi
import matplotlib.pyplot as plt
warnings.simplefilter("ignore", np.ComplexWarning)
def legendrePoly(a,n):
def integrand(t):
return ( ((a + exp(2j*pi*t))**2 - 1)/(2*exp(2j*pi*t)) )**n
return sp.integrate.quad(integrand,0,1)[0]
basisDim = 6
legendreBasis = [functools.partial(legendrePoly, n=i) for i in range(basisDim)]
integrand = [lambda x,i=i: exp(x) * legendreBasis[i](x) for i in range(basisDim)]
normalizingConst = [lambda x,i=i: legendreBasis[i](x)**2 for i in range(basisDim)]
basisCoeff = [sp.integrate.quad(integrand[i],-1,1)[0]
/sp.integrate.quad(normalizingConst[i],-1,1)[0] for i in range(basisDim)]
approxPoly = lambda x: sum(basisCoeff[i]*legendreBasis[i](x) for i in range(basisDim))
t = np.arange(-1, 1, 1e-3)
plt.plot(t,exp(t),'b')
plt.plot(t,approxPoly(t),'r')
plt.show()
I'm using Legendre polynomials as a basis for a polynomial approximation to the exponential function. I'm also using Cauchy's Integral Formula to evaluate them, rather than importing them directly from numpy.
Everything runs fine up to and including defining approxPoly, and approxPoly returns the expected value for any input I enter. But for some reason when I try to plot approxPoly(t), it returns the error: Supplied function does not return a valid float.
This error seems to suggest that when my functions in legendreBasis call to scipy.integrate.quad, that something goes wrong there, but if that was the case then approxPoly wouldn't work, yet if you evaluate it at 2000 points between -1 and 1 manually, and plot these points, everything works fine, but is that not exactly what plt.plot doing when it tries to graph my function?
Traceback:
Traceback (most recent call last):
File "/private/var/folders/mb/yyp8v3_95l538z3g7jsttq540000gn/T/Cleanup At Startup/Exercise-413072597.643.py", line 34, in <module>
plt.plot(t,approxPoly(t),'r')
File "/private/var/folders/mb/yyp8v3_95l538z3g7jsttq540000gn/T/Cleanup At Startup/Exercise-413072597.643.py", line 29, in <lambda>
approxPoly = lambda x: sum(basisCoeff[i]*legendreBasis[i](x) for i in range(basisDim))
File "/private/var/folders/mb/yyp8v3_95l538z3g7jsttq540000gn/T/Cleanup At Startup/Exercise-413072597.643.py", line 29, in <genexpr>
approxPoly = lambda x: sum(basisCoeff[i]*legendreBasis[i](x) for i in range(basisDim))
File "/private/var/folders/mb/yyp8v3_95l538z3g7jsttq540000gn/T/Cleanup At Startup/Exercise-413072597.643.py", line 18, in legendrePoly
return sp.integrate.quad(integrand,0,1)[0]
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 247, in quad
retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/integrate/quadpack.py", line 312, in _quad
return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
quadpack.error: Supplied function does not return a valid float.
logout
Upvotes: 1
Views: 1596
Reputation: 282198
scipy.integrate.quad
doesn't integrate functions that return arrays. That means when you try to call approxPoly(t)
, t
gets passed around until it ends up in this function:
def legendrePoly(a,n):
def integrand(t):
return ( ((a + exp(2j*pi*t))**2 - 1)/(2*exp(2j*pi*t)) )**n
return sp.integrate.quad(integrand,0,1)[0]
integrand
returns an array, and sp.integrate.quad
chokes. The same problem occurs a number of times in your code. Everything seems to be written assuming scalar arguments.
You can probably fix this by calling vectorize
on approxPoly
:
plt.plot(t,np.vectorize(approxPoly)(t),'r')
NumPy will then call approxPoly
on each element of t
separately.
Upvotes: 2