Python/Scipy集成数组

Question

我正在尝试编写一个程序来执行以下操作：

• 从数组中获取V值
• 将V值转换成关于E的积分
• 将积分结果输出到数组I
• 图一对五

这个方程看起来很糟糕，但是除了V.Here is the equation以外，所有的东西都是常数。这个等式不是很重要。在

这个问题我该怎么办？我的尝试（如下所示）不计算从文件中读取的每个V值的积分。在

from scipy import integrate #integrate.quad
from numpy import *
import pylab
import datetime
import time
import os
import math

# import V
fn = 'cooltemp.dat'
V = loadtxt(fn,unpack=True,usecols=[1])

# variables
del1, del2, R, E, fE, fEeV = 1,2,1,2,1,1
e = 1.602176565*10**-19

# eqn = dint(abc)
a = E/( math.sqrt( E**2 - del1**2 ) )
b = ( E+ e*V )/( math.sqrt( ( E + e*V )**2) - del2**2)
c = fE-fEeV
d = 1/(e*R) # integration constant
eqn = a*b*c

# integrate 
result = quad(lambda E: eqn,-inf,inf)

# current
I = result*d

# plot IV curve
pylab.plot(V,I,'-r')

## customise graph
pylab.legend(['degree '+str(n),'degree '+str(q),'data'])
pylab.axis([0,max(x),0,max(y)])
pylab.xlabel('voltage (V)')
pylab.ylabel('current (A)')
tc = datetime.datetime.fromtimestamp(os.path.getmtime(fn))
pylab.title('IV curve\n'+fn+'\n'+str(tc)+'\n'+str(datetime.datetime.now()))
pylab.grid(True)
pylab.show()

*更新尝试：

^{pr2}$

Answer 1

你有几个问题：

传递给quad的“函数”总是返回eqn，它只是一个预先计算好的数字。您需要定义一个适当的函数，它将E的给定值作为输入并返回被积函数。此函数还需要为V假定一个固定值。假设您提供的代码计算给定值V和E的正确数量（我没有检查，只是复制粘贴）：

# import V
fn = 'cooltemp.dat'
V = loadtxt(fn,unpack=True,usecols=[1])
# print V

@np.vectorize
def result(x):
    def integrand(E):
        del1, del2, R, fE, fEeV = 1.0,2.0,1.0,1.0,1.0
        e = 1.602176565*10**-19
        a = E/( math.sqrt( E**2 - del1**2 ) )
        b = ( E+ e*x )/( math.sqrt( ( E + e*x )**2) - del2**2)
        c = fE-fEeV
        d = 1/(e*R) # integration constant
        return a * b * c
    return quad(integrand, -inf, inf)

I = result(V)

总结一下：

• result(v)计算固定值v的全积分（E上）
• integrand(E)在固定的E（积分变量）和固定的V（它从函数外部获取值，这就是被积函数的定义嵌套在result定义中的原因）
• @np.vectorize技巧只是一个很好的方便函数，它允许您将V的数组传递到result。Numpy将为您循环这些值，并返回一个数组而不是标量

Answer 2

您应该使用np.vectorize将数组传递到表达式中，然后将数组取回。 例如，这将计算以下公式（如果您好奇，则计算移动距离…）：

import numpy as np
from scipy.integrate import quad

spl=299792458.0    #speed of light in m/s
Mpc=3.0856E22      # Mpc in m
pc=3.0856E16       # pc in m
def HubbleTime(H0):
    return 3.0856e17/(H0/100.0)
def HubbleDist(H0):
        """returns the Hubble Distance (in Mpc) for given H_0"""
    return spl*HubbleTime(H0)/Mpc
def Integrand(z, Om, OLam):
        """ This is the E(z) function from Hogg (2000)
        Integrand(z, Om, OLam)
    """
    return ( Om*(1+z)3 + OLam)(-0.5)
def CosmComDist(z, H0=70, Om=0.30, OLam=0.70):
    """Gives the comoving distance at redshift z
    CosmComDist(z, H0=70, Om=0.30, OLam=0.70)
    """
    CMD=HubbleDist(H0)*quad(Integrand, 0, z, args=(Om, OLam))[0]
    return CMD
CosmComDist=np.vectorize(CosmComDist)
redshifts = np.linspace(0,1,100)
distances = CosmComDist(redshifts)