找到数值积分的根

Question

我试图用Python重现这个Mathematica程序:

它找到了数值积分的根,并形成了这些值的图.但是,我不能让我试图跑.

目前的尝试:

来自scipy.integrate import quad from scipy import integration 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

正如评论中已经提到的（@Hristo Iliev 和 @Pavel Annosov），quad 返回一个 tuple of stuff。如果您假设积分没有问题，就像您在 Mathematica 中所做的那样（尽管这不是一个好主意），那么在这个元组中您只需要第一个元素，它应该是积分结果。

但这只会给你一个数字，而不是 的函数T。要获得后者，您需要自己定义相应的函数，就像您在 Mathematica 中使用您的函数一样\Delta[T_]:=...

以下是一些可以帮助您入门的信息：

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

请注意，您需要使用args=(T,)语法将参数发送T到正在集成的函数：quad对函数的第一个参数进行集成，并且需要其他参数才能计算f(E, T)。

现在您可以将其提供gap(T)给fsolve，它也需要一个函数callable（更准确地说是a ）。

在稍微更一般的层面上，您不应该使用玻尔兹曼常数、hbar 等的数值（甚至 Mathematica 也抱怨！）。相反，你应该做什么，你应该以无量纲形式编写公式：以能量单位测量温度（因此 k_B = 1）等，在积分中进行适当的替换，以便处理无量纲参数的无量纲函数 --- 和然后让计算机处理这些。