Удельная теплоотдача с использованием модели Дебая + Эйнштейна в matplotlib - PullRequest
Новая
       19

Удельная теплоотдача с использованием модели Дебая + Эйнштейна в matplotlib

1 голос
/ 23 мая 2019

Я пытаюсь найти соответствие конкретным данным о теплоте, используя гамма T + m Debye_model + (1-m) * модель Эйнштейна, как указано ниже.

Cel + ph (T) = γ T + [αCDebye (T) + (1 - α) CEinstein (T)]

, где модели Дебая и Эйнштейна задаются уравнением3 и 4 во вложении.

Я попробовал следующий код в блокноте jupyter, следуя некоторым примерам в Интернете, но я не представляю, как можно объединить эти функции вместе, чтобы выполнить подбор.

Данные связаны https://www.dropbox.com/s/u0r2m3zwl8w77at/HC_ScPtBi.dat?dl=0 Столбец 1 - это температура, а столбец 3 - это данные Y, представляющие интерес.Модель в https://www.dropbox.com/s/9452fq7eydajr5o/Debye.pdf?dl=0

Код в https://www.dropbox.com/s/hk9b1t0agvt36zn/Untitled2.ipynb?dl=0

from matplotlib import pyplot 
import numpy as np
from scipy import integrate
from scipy.optimize import curve_fit
from scipy.integrate import quad

data=np.genfromtxt('HC_ScPtBi.dat', skip_header=1)
R=8.314
n=3
M=1
T=data[10:290,0]
c=data[10:290,2]
def plot_data():
    pyplot.scatter(T, c)
    pyplot.xlabel('$T [K]$')
    pyplot.ylabel('$C$')

plot_data()
def c_einstein(T, T_E):
    x = T_E / T
    return 3 *n*R*x**2 * np.exp(x) / (np.exp(x) - 1)**2

popt0, pcov0 = curve_fit(c_einstein, T, c, 250)
T_E = popt0[0]
delta_T_E = np.sqrt(pcov0[0, 0])
print(f"T_E = {T_E:.5} ± {delta_T_E:.3} K")
print(popt0)

plot_data()
#temps = np.linspace(10, T[-1], 100)
pyplot.plot(T, c_einstein(T, *popt0));

def integrand(y):
    return y**4 * np.exp(y) / (np.exp(y) - 1)**2

@np.vectorize
def c_debye(T, T_D):
    x = T / T_D
    return 9 *n*R*x**3 * quad(integrand, 0, 1/x)[0]

popt1, pcov1 = curve_fit(c_debye, T, c, 150)
T_D = popt1[0]
delta_T_D = np.sqrt(pcov1[0, 0])
print(f"T_D = {T_D:.5} ± {delta_T_D:.3} K")
print(popt1)
plot_data()
pyplot.plot(T, c_einstein(T, *popt0), label='Einstein')
pyplot.plot(T, c_debye(T,  *popt1), label='Debye')
pyplot.legend();

1 Ответ

0 голосов
/ 23 мая 2019

Если это может быть полезным, я получил превосходное соответствие модифицированному уравнению пика Вейбулла, где R-квадрат = 0,99999 и RMSE = 0,06878.

enter image description here

def Peak_WeibullPeak_Modified_model(x): # from zunzun.com
    a = 6.4654735487019195E+01
    b = 3.4517137038577323E+02
    c = -1.5940608784806631E+00
    d = 2.7331145870203617E+00

    return = a * numpy.exp(-0.5 * numpy.power(numpy.log(x/b) / c, d))
...