Numerical integration for modelling curve for superconductors (Python)

Question

I am a physicist who is trying to model the current-voltage characteristics of a superconductor-superconductor junction.

The equation for this model is:

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

Current ($I$ or I in the code) values are calculated by evaluating this integral for given voltages ($V$, or v in the code).

I have attempted this in Python. The code is shown below.

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

However, I receive OverflowError: math range error. Does anybody have any ideas how this can be overcome? Apologies for the 10**n and long integrals. The code runs when the exponentials are removed (returns 0), and herein lies the problem.

Any ideas how this can be modelled in Python, or any other language?

Answer 1

First, it's always good to further debug the issue and seeing exactly where (what term, for what parameters) the overflow is coming from. This was not clear to me from the question.

Once you know precisely what the problem is, you can better diagnose the issue. For example, lets take the overflow issue. This is somewhere above 1e300 if I remember correctly.

1. You should be asking yourself 'should I really be hitting a range where this variable is getting this high'?
2. If the answer is 'yes, but it is divided out or canceled by another term', you need to rewrite the equations in a way that combines those terms. That is, if you are calculating $(x +y)/z$ with $x,y,z$ huge numbers try defining $x'=x/z$ and $y'=y/z$ and integrate those variables. The goals here should be to write your equations so that you are always be comparing terms of similar magnitudes.
3. If you are convinced you need very large numbers for the problem, try converting your equations to a form in which the fundamental variables are in logspace.

Note that you're explicitly asking the code to do an improper integral (-inf,inf), so you need to make sure the integrand you're putting into it is well behaved for numerics along that whole domain! It may help to try truncating the integral at the smallest value you can, lets call it $E_{max}$, where you know physically that you don't expect any contribution to the integral from $E_{max}$ to inf.

Good luck!