Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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)

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 = array(I)

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 = array(I2)

pl.xlabel(r'Voltage ($V$)')
pl.ylabel(r'Current ($A$)')
pl.title('Theoretical I(V) curve')

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?

share|improve this question
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). – GertVdE Feb 27 '13 at 13:44

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!

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.