Find root of numerical integration

Question

I am trying to reproduce this Mathematica program in Python:

It finds the roots of an numerical integration, and forms a plot of these values. However, I cannot get my attempt to run.

Current attempt:

from scipy.integrate import quad from scipy import integrate from scipy.optimize import fsolve import pylab as pl import numpy as np

# Variables.
boltzmann_const = 1.38e-23
planck_const = 6.62e-34
hbar = planck_const / ( 2 * np.pi )
transition_temp = 9.2
gap_energy_at_zero_kelvin = 3.528 / ( 2 * transition_temp * boltzmann_const )
debye_freq = ( 296 * boltzmann_const ) / hbar

# For subtracting from root_of_integral
a_const = np.log( ( 1.13 * hbar * debye_freq ) / ( boltzmann_const * transition_temp) )
# For simplifying function f.
b_const = ( hbar * debye_freq ) / ( 2 * boltzmann_const)


def f( coherence_length, temp ):
    # Defines the equation whose integral will have its roots found. Epsilon = coherence length. Delta = Gap energy.

    squareRoot = np.sqrt( coherence_length*coherence_length + gap_energy*gap_energy )
    return np.tanh( ( ( b_const / temp ) * squareRoot ) / squareRoot )


def integrate( coherence_length, temp ):
    # Integrates equation f with respect to E, between 0 and 1. 

    return integrate.quad( f, 0, 1, args = ( temp, ) )[0]


def root_of_integral( temp ):
    # Finds the roots of the integral with a guess of 0.01.

   return fsolve( integrate, 0.01, args = ( temp, ) )


def gap_energy_values( temp ):
    # Subtracts a_const from each root found, to obtain the gap_energy_values.

    return root_of_integral( temp ) - a_const

Answer 1

Just as has been mentioned in comments already (by @Hristo Iliev and @Pavel Annosov), quad returns a tuple of stuff. If you're assuming the integration has no problems, as you seem to be doing in Mathematica (which is not that good of an idea though), out of this tuple you only need the first element, which is supposed to be the integration result.

But this will only give you a single number, not a function of T. To get the latter, you'll need to define the corresponding function yourself, just like you did in Mathematica with your \Delta[T_]:=...

Here are a few bits to get you started:

def f(E, T):
    """To be integrated over T."""
    temp = np.sqrt(E * E + T * T)
    return np.tanh(1477.92 * temp) / temp


def gap(T):
    """Integration result: \Delta(T)"""
    return quad(f, 0, 1, args=(T, ))[0]  #NB: [0] select the 1st element of a tuple

Notice that you need to use the args=(T,) syntax to send the T parameter to the function being integrated: quad integrates over the first argument of the function, and it needs other arguments to be able to evaluate f(E, T).

Now you can feed this gap(T) to fsolve, which also expects a function (a callable, to be more precise).

On a slightly more general level, you should not be using numerical values of the Boltzmann constant, hbar and the like (even Mathematica complains!). What you should do instead, you should write your formulas in the dimensionless form: measure temperature in energy units (hence k_B = 1) etc, do appropriate substitutions in the integrals, so that you deal with dimensionless functions of dimensionless arguments --- and then make a computer handle those.

Answer 2

This line:

integral = (integrate.quad(lambda E: np.tanh(1477.92*np.sqrt(E**2+x**2))/np.sqrt(E**2+x**2), 0, 1)

has unbalanced parentheses:

integral = integrate.quad(lambda E: np.tanh(1477.92*np.sqrt(E**2+x**2))/np.sqrt(E**2+x**2), 0, 1)

Would be much easier to see if you broke it up, e.g.

x_values = arange(0.01, 0.1, 0.0001)

delta = []
for x in x_values:
    def fun(E):
        distance = np.sqrt(E * E + x * x)
        return np.tanh(1477.92 * distance) / distance

    integral = integrate.quad(fun, 0, 1)
    delta_val = fsolve(integral, 1e-23) - a
    delta.append(delta_val)

pl.plot(delta, x_values)