# Numerical integration for modelling curve for superconductors (Python)

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?

• Please, as jeffdk mentions in his answer below: check your values -> E goes from -inf to +inf while you are subtracting a value g = $O(1e-45)$ from the square of E! And for debugging purposes, it might be a good idea to define a function F defining your integrand (instead of using the lambda form) so you can split it up in its factors/terms/parts to identify the culprit (although in this case, all your factors are in trouble). Commented Feb 27, 2013 at 13:44

## 1 Answer

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!