寻找数值积分的根

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

- 8765674

---

解释器提供的错误信息可能很有帮助？ - Chris Seymour

已根据要求进行编辑，谢谢！ - 8765674

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)

- Pavel Anossov

---

感谢您迄今为止的回答！这看起来很棒，但我收到了错误：http://pastie.org/private/frbdw2xsurchme6s5czk1w - 8765674

这可能应该是一个单独的问题，但是integrate.quad返回一堆东西的元组，而fsolve期望一个函数。由于我不擅长数学，也无法解析您的Mathematica程序，所以我不知道如何修复它。 - Pavel Anossov

在您的情况下，integral 可能是一对数字：(result, absoluteError)。您是否期望 integrate 返回一个函数？ - Pavel Anossov

你需要定义一个函数来评估积分，然后将该函数传递给根查找器 fsolve。 - Hristo Iliev

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中，它也期望一个函数（更确切地说是一个callable）。

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

- ev-br

---

是的，Hristo，“自然单位”可能是更好的措辞，也可能不是。我个人更喜欢将其视为在甚至触摸键盘之前对积分等进行适当替换 - 但这是个人偏好。看起来我们都同意，在代码中像kB = 1.23e-23这样的常数是不正确的。至于稳定性，我会先等待OP解决单位问题，然后我们就能看到是否存在数值稳定性的真正问题。 - ev-br

我同意通常使用无量纲量可以简化问题并提高数值稳定性。但在这种特殊情况下不是这样的。 - Hristo Iliev

你知道那个神奇的数字1477在适当的单位下会变成什么吗？我不知道。一旦集成修复好了，就可以使用fsolve，它的起始点是1e-23（我猜）。 - ev-br

神奇的数字1477.92是戴拜截止能量和~0.1K热能之比的一半。无论使用适当或不适当的单位，它都不会改变，除非积分形式发生变化。由于我不打算将其转化为有关铌的BCS理论的讲座，并且考虑到numpy.tanh()在大参数下可以很好地评估，让我们等待并看看OP会想出什么。迄今为止，他对他的“Python超导性”问题都没有提供任何合理的反馈... - Hristo Iliev

其实不需要赏金。好吧，错误信息非常明确：您需要在调用作用域中定义“T”的值。该语言不使用连续变量，它需要数字--- 您需要定义“T”的数字。 - ev-br