python scipy curve-fitting scipy-optimize quad

Curve-fitting a non linear equation

I'm trying to fit my thermal conductivity into the Debye-Callaway equation. However, one of my parameters is coming back negative. I've tried different initial guesses. So I'm attaching a code with thermal conductivity data from literature and the values that they got to understand where my model is going wrong. What concerns me is it's not very sensitive to initial parameters. I'd really appreciate some help regarding this.

Thank you so much!

The parameters that the paper reported were

A=6.73E-43

B=5.38E-18

The parameters I got:

Fitted parameters: 1.417070339751493e-44, B = 1.075485244260637e-15.

The functions are defined exactly the same as in paper:

Curve fitting

EDIT: I've added a second piece of code that still returns the second factor negatively. I'm not sure if it's because the factor is getting stuck at a local minimum somewhere?

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

# Constants
kb = 1.380649e-23  # Boltzmann constant (J/K)
hbar = 1.0545718e-34  # Reduced Planck's constant (J*s)
theta_D = 280  # Debye temperature (K)
v = 2300  # Average phonon group velocity (m/s) -


# Define the inverse relaxation times
def tau_PD_inv(omega, A):
    return A * omega**4

def tau_U_inv(omega, B, T):
    return B * omega**2 * np.exp(-theta_D / (3 * T))

def tau_GB_inv(omega):
    D = 0.7e-6  # Grain boundary characteristic length (m)
    return v / D

def tau_total_inv(omega, A, B, T):
    return tau_PD_inv(omega, A) + tau_U_inv(omega, B, T) + tau_GB_inv(omega) 

# Debye Callaway model for thermal conductivity
def kappa_lattice(T, A, B):
    def integrand(x):
        omega = x * kb * T / hbar
        tau_total = 1 / tau_total_inv(omega, A, B, T)
        return x**4 * np.exp(x) / (np.exp(x) - 1)**2 * tau_total

    integral, _ = quad(integrand, 0, theta_D / T)
    prefactor = kb**4 * T**3 / (2 * np.pi**2 * v * hbar**3)
    return prefactor * integral

# Vectorize the model for curve fitting
def kappa_lattice_vectorized(T, A, B ):
    return np.array([kappa_lattice(Ti, A, B ) for Ti in T])

# Example thermal conductivity and temperature data for Ge WQ
T_data = np.array([3.0748730372823156, 3.1962737114707203, 3.4441443065638015, 3.785528026732028,
    4.043999986888003, 4.524033859711473, 4.902329259251676, 5.270E+00,
    5.902534836759174, 6.721645553370004, 7.690683028568884, 9.031E+00,
    1.103E+01, 12.811754075061538, 16.07373842258417, 19.67675447459597,
    23.951202854256202, 29.710517518825778, 3.684E+01, 45.030707655028536,
    56.92469139605285, 6.929E+01, 85.62796640159988, 105.41803279520694,
    129.5096325800293, 159.77670578646953, 200.9314386671018, 240.34130109997898,
    291.63127697234717])  # Temperature in K
kappa_data = np.array([ 0.5021456742832973, 0.5907837911587943, 0.7331309897165834, 0.9570228747931655,
    1.1281866396926783, 1.4904675008427906, 1.8682223847710366, 2.1856825810901013,
    2.8515011626042894, 4.057791636291224, 5.328111104734368, 7.645627510566391,
    11.045409412114584, 14.277055520040301, 19.54065663430087, 24.738891215975602,
    28.776422188514683, 32.562485211039274, 33.37572711709385, 32.455084248250884,
    28.56802589740088, 24.937473590597776, 20.887192910896978, 17.37073097716312,
    14.783541083880793, 12.26010669954394, 9.991973525044475, 8.531955488387288,
    7.347906481962733])  # Thermal conductivity in W/mK

# Initial guess for parameters A, B, and C
initial_guess = [4.75E-42, 1.054E-17]

# Curve fitting
popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess, maxfev=8_000)
#popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess, bounds=(0, np.inf))

# Extract fitted parameters
A_fit, B_fit = popt

# Print the fitted parameters
print(f"Fitted parameters: A = {A_fit}, B = {B_fit}")

# Generate curve fit data
T_fit = np.linspace(min(T_data), max(T_data), 100)
kappa_fit = kappa_lattice_vectorized(T_fit, *popt)


# Plot the original data, fitted data, and extrapolated data
plt.figure(figsize=(8, 6))
plt.plot(T_data, kappa_data, 'o', label='Experimental data')
plt.plot(T_fit, kappa_fit, '-', label='Fitted model')
#plt.plot(T_extrap, kappa_extrap, '--', label='Extrapolated (400-600K)')
plt.xlabel('Temperature (K)')
plt.ylabel('Thermal Conductivity (W/mK)')
plt.legend()
plt.title('Thermal Conductivity Fitting ')
plt.show()

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

# Constants
kb = 1.380649e-23  # Boltzmann constant (J/K)
hbar = 1.0545718e-34  # Reduced Planck's constant (J*s)
theta_D = 483.0  # Debye temperature (K)
v = 4618.0  # Average phonon group velocity (m/s) -


# Define the inverse relaxation times
def tau_PD_inv(omega, A):
    return A * omega**4

def tau_U_inv(omega, B, T):
    return B * omega**2 * T * np.exp(-theta_D / (3 * T))

def tau_GB_inv(omega):
    D = 1e-6  # Grain boundary characteristic length (m)
    return v / D

def tau_new_inv(omega, C):
    return C * omega**2

def tau_total_inv(omega, A, B, C, T):
    return tau_PD_inv(omega, A) + tau_U_inv(omega, B, T) + tau_GB_inv(omega) + tau_new_inv(omega, C)

# Debye Callaway model for thermal conductivity
def kappa_lattice(T, A, B, C):
    def integrand(x):
        omega = x * kb * T / hbar
        tau_total = 1 / tau_total_inv(omega, A, B, C, T)
        return x**4 * np.exp(x) / (np.exp(x) - 1)**2 * tau_total

    integral, _ = quad(integrand, 0, theta_D / T)
    prefactor = kb**4 * T**3 / (2 * np.pi**2 * v * hbar**3)
    return prefactor * integral

# Vectorize the model for curve fitting
def kappa_lattice_vectorized(T, A, B, C):
    return np.array([kappa_lattice(Ti, A, B, C) for Ti in T])

# Example thermal conductivity and temperature data for Ge WQ
T_data = np.array([31.21198032, 31.63863233, 33.88415583, 36.02698918, 38.39872854,
    39.56437842, 43.09065509, 44.50006634, 46.62080745, 43.20705304,
    47.77194375, 51.1374173, 53.33894737, 55.68342411, 57.7839973,
    60.44084481, 63.31137245, 66.34735373, 69.16983511, 74.09615366,
    77.63998236, 82.66616204, 87.95740211, 93.19748533, 99.59297955,
    105.6903923, 113.1334233, 120.6776065, 127.0279836, 134.1693093,
    141.3911575, 148.1783285, 155.5512459, 162.631427, 169.7085111,
    177.378102, 183.5898785, 190.9367188, 197.6917958, 205.7338436,
    212.8065032, 220.0698325, 226.7609477, 234.0571384, 241.1287763,
    247.8786187, 255.2722521, 262.3425521, 269.4163916, 276.4872816,
    283.5510925, 297.1613034, 304.7919793, 311.8668132, 318.9352276,
    326.0927208, 333.0893085, 340.164332, 347.2375209, 354.9542919,
    362.0246034, 369.0955138, 375.8247965, 383.5799572])  # Temperature in K
kappa_data = np.array([ 1.085444844, 1.13592882, 1.176429149, 1.217086899,
    1.258678341, 1.305290878, 1.357081371, 1.398857871, 1.26084173,
    1.436622858, 1.483042548, 1.528031021, 1.561551609, 1.598710119,
    1.632418606, 1.665656346, 1.703762245, 1.742520105, 1.774186181,
    1.820109169, 1.855169951, 1.893680154, 1.932568555, 1.962172074,
    1.999232142, 2.029595098, 2.0490363, 2.067859052, 2.092450212,
    2.119053495, 2.131413546, 2.145488003, 2.158842932, 2.169084281,
    2.169487956, 2.184424479, 2.196747698, 2.20620033, 2.217744549,
    2.223537925, 2.236483003, 2.258522251, 2.273748326, 2.278514432,
    2.28270469, 2.288247862, 2.291668986, 2.298648488, 2.302662675,
    2.299560106, 2.311410793, 2.339970141, 2.347949378, 2.349474819,
    2.354437976, 2.369869851, 2.378039597, 2.384365082, 2.39183344,
    2.395266095, 2.399300875, 2.40766296, 2.409888354, 2.442162261])  # Thermal conductivity in W/mK

# Initial guess for parameters A, B, and C
initial_guess = [4.75E-43, 3.6979E-19, 1e-16]

# Curve fitting
popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess,maxfev=8_000)
#popt, pcov = curve_fit(kappa_lattice_vectorized, T_data, kappa_data, p0=initial_guess, bounds=(0, np.inf))

# Extract fitted parameters
A_fit, B_fit, C_fit = popt

# Print the fitted parameters
print(f"Fitted parameters: A = {A_fit}, B = {B_fit}, C = {C_fit}")

# Generate curve fit data
T_fit = np.linspace(min(T_data), max(T_data), 100)
kappa_fit = kappa_lattice_vectorized(T_fit, *popt)


# Plot the original data, fitted data, and extrapolated data
plt.figure(figsize=(8, 6))
plt.plot(T_data, kappa_data, 'o', label='Experimental data')
plt.plot(T_fit, kappa_fit, '-', label='Fitted model')
#plt.plot(T_extrap, kappa_extrap, '--', label='Extrapolated (400-600K)')
plt.xlabel('Temperature (K)')
plt.ylabel('Thermal Conductivity (W/mK)')
plt.legend()
plt.title('Thermal Conductivity Fitting ')
plt.show()

Solution

Your tau_U_inv function forgot to take into account T.

The corrected version should be this:

def tau_U_inv(omega, B, T):
    return B * omega**2 * T * np.exp(-theta_D / (3 * T))

When I make this change the parameters show

Fitted parameters: A = 1.935907419814169e-43, B = 9.770738747258265e-18

which appears to be much closer to the original.