6

私はこのMathematicaプログラムをPythonで再現しようとしています：

Mathematicaプログラム

数値積分の根を見つけ、これらの値のプロットを形成します。しかし、実行しようとはできません。

現在の試み：

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

2 に答える 2

2

この行：

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)
于 2013-02-22T15:30:36.203 に答える
2

すでにコメントで言及されているように（@HristoIlievと@PavelAnnosovによって）、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)をにフィードできます。これには、関数（より正確にはa）も必要です。fsolvecallable

もう少し一般的なレベルでは、ボルツマン定数やhbarなどの数値を使用するべきではありません（Mathematicaでさえ文句を言います！）。代わりに、式を無次元形式で記述する必要があります。エネルギー単位で温度を測定し（したがって、k_B = 1）、積分で適切な置換を行い、無次元引数の無次元関数を処理します---および次に、コンピューターにそれらを処理させます。

于 2013-02-22T18:59:44.347 に答える