python - Python/Scipy 統合配列

Question

次のことを行うプログラムを作成しようとしています。

配列から V の値を取ります
V 値を E に関する積分に渡します。
積分結果を配列に出力 I
I を V に対してプロットする

方程式はかなり厄介に見えますが、V 以外はすべて定数です。方程式はあまり重要ではありません。

この問題はどうすればいいですか？私の試み（以下に示すように）は、ファイルから読み取った 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()

*更新された試み：

from scipy import integrate
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])
# print V

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

I=[]
for n in range(len(V)):

    constant = 1/(e*R) # integration constant
    eqn = (E/( math.sqrt( E**2 - del1**2 ) ))*(( E + e*V[n] )/( math.sqrt( ( E + e*V[n] )**2) - del2**2))*(fE-fEeV)

    # integrate 
    result,error = integrate.quad(lambda E: eqn,-inf,inf)
    print result
    # current
    I.append(result*constant)

I = array(I)

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

score 4 · Accepted Answer

いくつかの問題があります。

渡す「関数」は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 (関数の外部から値を取得します。これが、被積分関数の定義が結果の定義内にネストされている理由です)
トリックは、@np.vectorizeV の配列をに渡すことができる便利な関数ですresult。Numpy は値をループして、スカラーではなく配列を返します。

score 1 · Accepted Answer

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)

python - Python/Scipy 統合配列

2 に答える 2

Related

Reference