如果我有个函数 f(x) = a*x^2+b*x+c      (a, b,c 为常数)

x = range(0,100)

如何求出 y(x) = 对alpha的积分 [  f(x)*cos(alpha) ], 其中alpha = [0, pi]

具体定义如下：

import numpy

from numpy import exp

def function1(x,w,xc):

    M_LN2 = 0.69314718055994530942

    yh = numpy.sqrt(M_LN2)*w[1]/w[0]

    xh = 2*numpy.sqrt(M_LN2)*(x-xc)/w[0]

    A =  2 * numpy.sqrt(M_LN2/numpy.pi)/w[0]

    y = exp(-A * function2(xh,yh))    

    return y

def function2(x,y):    

    T = [0.314240376254359,0.947788391240164,1.597682635152605,2.27950708050106,3.020637025120890,3.889724897869782]

    alpha = [-1.393236997981977,-0.231152406188676,0.155351465642094,-6.21836623696556e-3,-9.190829861057113e-5,6.275259577497896e-7]

    beta = [1.011728045548831,-0.751971469674635,1.255772699323164e-2,1.0022008145159e-2,-2.42068134815573e-4,5.008480613664573e-7]

    sum = numpy.zeros(x.shape,'double')

    y0 = 1.50

    yp = y + y0 

    yp2 = yp * yp

    if (y > 0.85):

        for k in range(0,5):

            xp = x+ T[k]

            xm = x - T[k]

            sum += ((alpha[k] * xm + beta[k] * yp) / (xm * xm + yp2) + (beta[k] * yp - alpha[k] * xp) / (xp * xp + yp2))

    indx =(numpy.absolute(x) < (18.1 * y + 1.65)).nonzero()

    #print indx

    sum[indx]  = 0

    for k in range(0,5):

        xp = x[indx] + T[k]

        xm = x[indx] - T[k]

        sum[indx] += ((alpha[k] * xm + beta[k] * yp) / (xm * xm + yp2) + (beta[k] * yp - alpha[k] * xp) / (xp * xp + yp2))

    indx2 = (numpy.absolute(x) >= (18.1 * y + 1.65)).nonzero()

    sum[indx2]=0

   # print indx2

    for k in range(0,5):

        xp = x[indx2] + T[k];

        xp2 = xp * xp;

        xm = x[indx2] - T[k];

        xm2 = xm * xm;

        sum[indx2] += (((beta[k] * (xm2 - y0 * yp) - alpha[k] * xm * (yp + y0))/ ((xm2 + yp2) * (xm2 + y0 * y0)))+ ((beta[k] * (xp2 - y0 * yp) + alpha[k] * xp * (yp + y0))/ ((xp2 + yp2) * (xp2 + y0 * y0))))

    indx3 = ((numpy.absolute(x) < 100.) & (numpy.absolute(x) >= (18.1 * y + 1.65))).nonzero()

    sum[indx3] =y*sum[indx3]+numpy.exp(-pow(x[indx3],2))

    indx4 = ((numpy.absolute(x) >= 100.) & (numpy.absolute(x) >= (18.1 * y + 1.65))).nonzero()

    sum[indx4] = y*sum[indx4]

    return sum

if __name__=="__main__":

    import pylab

    from numpy import cos, pi 

    from scipy import integrate

    from scipy.integrate import quad, dblquad

    x = numpy.arange(0,2048)

    y = function1(x,(6,100),1000) 

    # 先对 t = [0, pi] 积分后再计算出每个 x 对应的值。 

    #jifen = integrate.quad(y(x)*cos(t), 0, pi) # 错误的计算方法， 如何改正？？？

     dblquad(lambda t, x:   function1 *cos(t), 0,  pi, lambda x: 1, lambda x: 2048)

    pylab.figure()

    pylab.plot(x,y)

    pylab.plot(x,jifen)

    pylab.show()

TypeError: unsupported operand type(s) for *: 'function' and 'numpy.float64'

这个积分直接可以用 scipy.integrate.quad 做吧。

def y(x):
    return quad(lambda alpha: f(x)*cos(alpha), 0, pi)