Численное интегрирование для моделирования кривой для сверхпроводников (Python)

Я физик, который пытается смоделировать вольт-амперные характеристики соединения сверхпроводник-сверхпроводник.

Уравнение для этой модели:

\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}$

Текущий ( $I$ или Iв коде) значения рассчитываются путем оценки этого интеграла для заданных напряжений ( $V$ или vв коде).

Я пытался это в Python. Код показан ниже.

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()

Однако я получаю OverflowError: math range error. У кого-нибудь есть идеи, как это можно преодолеть? Извинения за 10**nи длинные интегралы. Код выполняется, когда экспоненты удалены (возвращает 0), и в этом заключается проблема.

Есть идеи, как это можно смоделировать в Python или любом другом языке?

python numpy integral-equations запрос
источник

Пожалуйста, как упоминает Джеффдк в своем ответе ниже: проверьте свои значения -> E переходит от -inf к + inf, пока вы вычитаете значение g =

O (1 e - 45)

$O(1e-45)$ с площади Е! А для целей отладки было бы неплохо определить функцию F, определяющую вашу подынтегральную функцию (вместо того, чтобы использовать лямбда-форму), чтобы вы могли разделить ее на ее факторы / термины / части, чтобы определить виновника (хотя в этом случае все ваши факторы в беде).

GertVdE

Ответы:

Во-первых, всегда полезно продолжить отладку проблемы и точно выяснить, откуда (какой срок, для каких параметров) происходит переполнение. Это было не ясно для меня из вопроса.

Как только вы точно знаете, в чем проблема, вы можете лучше диагностировать проблему. Например, давайте возьмем проблему переполнения. Это где-то выше 1е300, если я правильно помню.

Вы должны спросить себя: «Должен ли я действительно достигнуть диапазона, в котором эта переменная становится такой высокой»?
Если ответ «да, но он разделен или отменен другим термином», необходимо переписать уравнения таким образом, чтобы объединить эти термины. То есть, если вы рассчитываете $(x +y)/z$ с $x,y,z$ огромные числа пытаются определить $x'=x/z$ а также $y'=y/z$ и интегрировать эти переменные. Цели здесь должны состоять в том, чтобы написать свои уравнения так, чтобы вы всегда сравнивали схожие величины.
Если вы уверены, что вам нужны очень большие числа для задачи, попробуйте преобразовать ваши уравнения в форму, в которой фундаментальные переменные находятся в лог-пространстве.

Обратите внимание, что вы явно просите код сделать неправильный интеграл (-inf, inf), поэтому вам нужно убедиться, что интегрируемый в него код корректно ведет себя для чисел во всем этом домене! Это может помочь попытаться усечь интеграл до наименьшего возможного значения, давайте назовем его $E_{max}$ где вы физически знаете, что не ожидаете какого-либо вклада в $E_{max}$ до инф.

Удачи!

jeffdk
источник