# 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:
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). Feb 27 '13 at 13:44

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.
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.