6

モデルⅣ。

方法: E の関数として積分を実行し、使用される各電圧値の電流を出力します。これは、v_values の配列に対して繰り返されます。方程式は以下で見つけることができます。

ここに画像の説明を入力

この方程式の範囲は から-infまでinfですが、極を避けるために (E+eV)^2-\Delta^2>0 および E^2-\Delta^2>0 となるように制限する必要があります。(\Delta_1 = \Delta_2)。したがって、現在 2 つの積分があり、範囲は-infから-gap-e*v およびgapまでinfです。

math range errorただし、上記の制限を使用して厄介な E 値を除外したと思いますが、a を返し続けています。エラーの貼り付け: http://pastie.org/private/o3ugxtxai8zbktyxtxuvg

この質問のあいまいさをお詫びします。しかし、誰でも明らかな間違いやコードの誤用を見つけることができますか?

私の試み:

from scipy import integrate
from numpy import *
import scipy as sp
import pylab as pl
import numpy as np
import math

e = 1.60217646*10**(-19)
r = 3000
gap = 400*10**(-6)*e
g = (gap)**2
t = 0.02
k = 1.3806503*10**(-23)
kt = k*t

v_values = np.arange(0,0.001,0.0001)

I=[]
for v in v_values:
    val, err = integrate.quad(lambda E:(1/(e*r))*(abs(E)/np.sqrt(abs(E**2-g)))*(abs(E+e*v)/(np.sqrt(abs((E+e*v)**2-g))))*((1/(1+math.exp((E+e*v)/kt)))-(1/(1+math.exp(E/k*t)))),-inf,(-gap-e*v)*0.9)
    I.append(val)
I = array(I)

I2=[]
for v in v_values:
    val2, err = integrate.quad(lambda E:(1/(e*r))*(abs(E)/np.sqrt(abs(E**2-g)))*(abs(E+e*v)/(np.sqrt(abs((E+e*v)**2-g))))*((1/(1+math.exp((E+e*v)/kt)))-(1/(1+math.exp(E/k*t)))),gap*0.9,inf)
    I2.append(val2)
I2 = array(I2)

I[np.isnan(I)] = 0
I[np.isnan(I2)] = 0

pl.plot(v_values,I,'-b',v_values,I2,'-b')
pl.show()
4

3 に答える 3

4

この質問は、計算科学のサイトに適しています。それでも、ここであなたが考えるいくつかのポイントがあります。

まず、積分範囲は と の交点(-oo, -eV-gap) U (-eV+gap, +oo)です(-oo, -gap) U (gap, +oo)。次の 2 つのケースが考えられます。

  • eV < 2*gap許容されるエネルギー値が にある場合(-oo, -eV-gap) U (gap, +oo)
  • の場合eV > 2*gap、許容されるエネルギー値は です(-oo, -eV-gap) U (-eV+gap, -gap) U (gap, +oo)

第二に、あなたは非常に低温の地域で働いています。0.02 Kにt等しい場合、ボルツマン係数の分母は 1.7 µeV であり、エネルギー ギャップは 400 µeV です。この場合、指数の値は正のエネルギーに対して巨大であり、Python で使用される倍精度浮動小数点数の限界をすぐに超えてしまいます。これは可能な最小の正のエネルギーであるため、より高いエネルギーでは物事は良くなりません. 負のエネルギーでは、値は常にゼロに非常に近くなります。この温度では、フェルミ ディラック分布のエッジが非常に鋭く、反射されたシータ関数に似ていることに注意してください。で、約 6.24E+100 にE = gapなります。またはexp(E/kT)の場合、範囲外になります。E/kT > 709.78E > 3.06*gap

しかし、その温度では、 (0.8 mV)[-eV, 0]の場合、与えられた温度のギャップ内に完全に収まる間隔の外側で、2 つのフェルミ関数の差が非常に急速にゼロになるため、そのようなエネルギーに行くことは意味がありません。V < (2*gap)/eそのため、バイアス電圧が 0.8 mV 未満の場合、電流はゼロに非常に近くなると予想されます。0.8 mV を超える場合、積分の主な値は の被積分関数から得られますが、(-eV+gap, -gap)ゼロ以外の値の一部は での特異点E = gap付近の領域から、一部は での特異点付近の領域から得られE = -eV-gapます。DoS の特異点を避けるべきではありません。そうしないと、I(V) 曲線に予期される不連続点 (垂直線) が得られません (画像は から取得)。ウィキペディア):

STJ電流電圧図

むしろ、各特異点の近傍で同等の近似式を導出し、代わりにそれらを統合する必要があります。

ご覧のとおり、被積分関数の値には多くの特殊なケースがあり、数値計算を行う際にはそれらすべてを考慮に入れる必要があります。それをしたくない場合は、Maple や Mathematica などの他の数学パッケージを使用する必要があります。これらには、はるかに洗練された数値積分ルーチンがあり、数式を直接処理できる場合があります。

これはあなたの質問に答えようとするものではなく、コメント フィールドに収まらない非常に長いコメントであることに注意してください。

于 2013-02-19T11:33:14.597 に答える
2

数学範囲エラーの理由は、指数が無限大になることです。とを取るv = 0.0009E = 5.18e-23、式exp((E + e*v) / kt)(Python 式で Hristo Liev が指摘したタイプミスを修正しました) はexp(709.984..)、倍精度数 (最大で 1E308) で表現できる範囲を超えています。

2 つの追加メモ:

  • 他の人が指摘したように、より小さな範囲の数値を提供する単位系を使用して、方程式を再スケーリングする必要があります。おそらく、原子単位は設定されるので可能な選択ですe = 1が、私はあなたの方程式をそれに変換しようとしませんでした. (おそらく、原子単位では時間単位が約 1/40 fs であるため、タイムステップは非常に大きくなります)。

  • 通常、浮動小数点数には指数表記を使用します:e = 1.60217E-19の代わりにe = 1.60217*10**(-19).

于 2013-02-18T21:43:44.473 に答える
1

最終的にこの問題に対処する最善の方法は、heaviside 関数を使用して、E変数が変数を超えないようにすることでし\Deltaた。

于 2013-03-01T13:41:14.307 に答える