寻找数值积分的根

Question

我正在尝试用Python重现这个Mathematica程序：

这个程序的功能是找到数值积分的根，并绘制这些值的图形。不过，我的尝试一直无法运行。

目前的尝试代码：

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

这一行：

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

有不匹配的括号：

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

如果你把它拆开来看，会更容易理解，比如：

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)

Answer 2

正如评论中提到的（@Hristo Iliev 和 @Pavel Annosov），quad 返回的是一个元组。如果你认为积分没有问题，就像你在 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，它也需要一个函数（更准确地说，是一个 可调用对象）。

从更一般的角度来看，你不应该使用玻尔兹曼常数、hbar 等的数值（连 Mathematica 都会抱怨！）。你应该做的是，把你的公式写成无量纲的形式：用能量单位来衡量温度（因此 k_B = 1）等等，在积分中做适当的替换，这样你就能处理无量纲的函数和无量纲的参数——然后让计算机来处理这些。