Całkowanie numeryczne krzywej modelowania dla nadprzewodników (Python)

Jestem fizykiem, który próbuje modelować charakterystykę prądowo-napięciową złącza nadprzewodnik-nadprzewodnik.

Równanie dla tego modelu to:

\begin{aligned} I (V) = \frac{1}{e R_{n - n}} \int_{- \infty}^{\infty} \frac{| E |}{[E^{2} - Δ_{1}^{2}]^{1 / 2}} \frac{| E + e V |}{[(E + e V)^{2} - Δ_{2}^{2}]^{1 / 2}} [f (E) - f (E + e V)] d E \end{aligned}

$\begin{align} I(V) = \frac{1}{eR_{\mathrm{n-n}}}\int_{-\infty}^{\infty}\frac{|E|}{[E^{2} - \Delta_{1}^{2}]^{1/2}}\frac{|E + eV|}{[(E + eV)^{2} - \Delta_{2}^{2}]^{1/2}}[f(E) - f(E + eV)]\,\mathrm{d}E \end{align}$

Wartości prądu ( lub w kodzie) oblicza się, oceniając tę całkę dla danych napięć ( lub w kodzie). $I$ I $V$ v

Próbowałem tego w Pythonie. Kod pokazano poniżej.

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

ec = 1.6021764*10**(-19)
r = 2500
gap = 200*10**(-6)*ec
g = (gap)**2
t = 0.04
k = 1.3806503*10**(-23)
kt = k*t

v_values = arange(0,0.001,0.00001)

I=[]
for v in v_values:
    result, error = integrate.quad (lambda E:(abs(E)/sqrt((E**2-g)))*(abs(E+ec*v)/(sqrt(((E+ec*v)**2-g))))*(math.exp(-E/kt)*(math.exp(-ec*v/kt)-1)),(-inf),(-gap*0.9-ec*v))
    I.append(result)
I = array(I)

I2=[]
for v in v_values:
    result2 = integrate.quad(lambda E:(abs(E)/sqrt((E**2-g)))*(abs(E+ec*v)/(sqrt(((E+ec*v)**2-g))))*(math.exp(-E/kt)*(math.exp(-ec*v/kt)-1)),(gap*0.9),(inf))
    I2.append(result2)
I2 = array(I2)

pl.plot(v_values,I,'-b',v_values,I2,'-r')
pl.xlabel(r'Voltage ($V$)')
pl.ylabel(r'Current ($A$)')
pl.title('Theoretical I(V) curve')
pl.grid(True)
pl.savefig('IVcurve.png')
pl.show()

Jednak otrzymuję OverflowError: math range error. Czy ktoś ma jakieś pomysły, jak to rozwiązać? Przepraszamy za 10**ncałki i długie. Kod jest uruchamiany, gdy wykładnicze są usuwane (zwraca 0), i na tym polega problem.

Wszelkie pomysły, w jaki sposób można to modelować w języku Python lub innym języku?

python numpy integral-equations

— pytanie
źródło

Proszę, jak wspomina jeffdk w swojej odpowiedzi poniżej: sprawdź swoje wartości -> E przechodzi od -inf do + inf podczas odejmowania wartości g = od kwadratu E! A do celów debugowania dobrym pomysłem może być zdefiniowanie funkcji F definiującej twój integrand (zamiast używania postaci lambda), abyś mógł podzielić go na czynniki / warunki / części w celu zidentyfikowania sprawcy (chociaż w tym przypadku wszystkie twoje czynniki mają kłopoty).

O (1 e - 45)

$O(1e-45)$

— GertVdE

Po pierwsze, zawsze dobrze jest debugować problem i zobaczyć dokładnie, skąd (jaki termin, dla jakich parametrów) pochodzi przepełnienie. To pytanie nie było dla mnie jasne.

Gdy dokładnie wiesz, na czym polega problem, możesz lepiej zdiagnozować problem. Weźmy na przykład problem z przepełnieniem. To jest gdzieś powyżej 1e300, jeśli dobrze pamiętam.

Powinieneś zadać sobie pytanie „czy naprawdę powinienem trafić w zakres, w którym ta zmienna osiąga tak wysoki poziom”?
Jeśli odpowiedź brzmi „tak, ale jest podzielona lub anulowana przez inny termin”, musisz przepisać równania w sposób, który łączy te terminy. To znaczy, jeśli obliczasz pomocą wielkich liczb spróbuj zdefiniować oraz i zintegruj te zmienne. Celem powinno być napisanie równań, abyś zawsze porównywał terminy o podobnych wielkościach. $(x +y)/z$ $x,y,z$ $x'=x/z$ $y'=y/z$
Jeśli jesteś przekonany, że potrzebujesz bardzo dużych liczb do rozwiązania problemu, spróbuj przekształcić swoje równania do postaci, w której podstawowe zmienne znajdują się w przestrzeni logów.

Zauważ, że wyraźnie prosisz kod o wykonanie niepoprawnej całki (-inf, inf), więc musisz upewnić się, że integra, w którą wpisujesz, jest dobrze zachowana dla liczb w całej tej domenie! Pomoże to spróbować skrócić całkę przy najmniejszej możliwej wartości, nazwijmy to , gdzie fizycznie wiesz, że nie oczekujesz żadnego wkładu w całkę od do inf. $E_{max}$ $E_{max}$

Powodzenia!

— jeffdk
źródło