Find the root of numerical integration

Question

I am trying to reproduce this Mathematica program in Python:

He finds the roots of numerical integration and forms a graph of these values. However, I cannot complete my attempt.

Current attempt:

from scipy.integrate import quad from using scipy import 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

As mentioned in the comments (@Hristo Iliev and @Pavel Annosov), quad returns a tuple of things . If you assume that integration does not have problems, as you seem to do in Mathematica (which is not so good for an idea), from this tuple you need only the first element, which should be the result of integration,

But this will give you only one number, not a T function. To get the latter, you will need to define the corresponding function yourself, as in Mathematica using \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 

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

Now you can send this gap(T) to fsolve , which also expects a function (a callable , more precisely).

On a slightly more general level, you should not use the numerical values ​​of the Boltzmann constant, hbar, etc. (even Mathematica complains!). Instead, you should write your formulas in dimensionless form: measure temperature in units of energy (therefore, k_B = 1), etc., make appropriate replacements in the integrals, so that you are dealing with dimensionless functions of dimensionless arguments --- and then get the computer to process them.

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) 

It would be much easier to see if you broke it, for example.

 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)