Cannot figure out how to fit a complicated custom equation.

I am having trouble implementing this equation into the curve fitting toolbox, both in the command and by utilizing the curve fitting GUI.

A little background that may help explain this:

I am wanting to fit experimental data of T(K) vs K_ph (W/m K) which the thermal conductivity of a material, to a model from a journal paper. I have done an okay job in python with scipy.optimize but have been trying MATLAB for a more accurate method.

I have this equation I need to fit

where

A, B, b, C_1, C_2, D, Omega_res1, Omega_res2, and L (ignore L in the exmaple code) are all coefficients I want solved for.

Omega is the independent variable in the integral (omega is solved for) while T is the range 1 to 100 in 1K steps. Or T can simply be my experimental Temperature data as is in my code.

My problems (in the GUI) are in the ways to set the values for T and use an integral in the custom equation box. As in, I am not sure how to tyep the equation in the way MATLAB accepts in the GUI for the toolbox.

My problems for the code is in the integration and lsqcurve fit it seems. I have a lot of warnings/errors such as

%{
Warning: Derivative finite-differencing step was artificially reduced to be within bound constraints. This may
adversely affect convergence. Increasing distance between bound constraints, in dimension 6, to be at least 2e-20
may improve results. 

Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.

<stopping criteria details>
Fitted Parameters:
D = 1e-40
B = 1e-20
A = 1e-15
b = 100
C1 = 1e-40
C2 = 1e-40
%}

I am not sure how to fix these, I have gone through other posts on similar matters and have changed the bounds and manually changed the optimoptions values to very low, but reached a point where it was just iterating nonstop.

Any help is appreciated!

% Data from 0T.txt
data = readmatrix('0T.txt', 'HeaderLines', 1);
T_data = data(:, 1); % First col: Temperature (K)
k_ph_data = data(:, 2); % Second col: \kappa_{xx} (W/mK)

% Constants
k_B = 1.380649e-23; % Boltzmann constant in J/K
hbar = 1.0545718e-34; % Reduced Planck constant in J·s
v_s = 4.817e3; % Sound velocity in m/s
theta_D = 235; % Debye temperature in K
L = 1e-6; % Length in m

% Define the integrand function with omega_res1 and omega_res2
integrand = @(omega, T, D, B, A, b, C1, C2, omega_res1, omega_res2) ...
    (omega.^4 .* exp(hbar*omega./(k_B*T)) ./ (exp(hbar*omega./(k_B*T)) - 1).^2) .* ...
    (v_s / L + D*omega.^4 + B*omega + A*T*omega.^3 .* exp(-theta_D / (b*T)) + ...
    C1 * omega.^4 ./ ((omega.^2 - omega_res1.^2).^2 + 1e-10) .* exp(-hbar*omega_res1 / (k_B * T)) ./ ...
    (1 + exp(-hbar*omega_res1 / (k_B * T))) + ...
    C2 * omega.^4 ./ ((omega.^2 - omega_res2.^2).^2 + 1e-10) .* exp(-hbar*omega_res2 / (k_B * T)) ./ ...
    (1 + exp(-hbar*omega_res2 / (k_B * T))));

% Define the k_ph function
k_ph_func = @(params, T) (k_B / (2 * pi^2 * v_s)) * (k_B * T / hbar)^3 .* ...
    integral(@(omega) integrand(omega, T, params(1), params(2), params(3), params(4), params(5), params(6), params(7), params(8)), 0, theta_D / T);

% Define the function to be minimized
fun = @(params, T) arrayfun(@(t) k_ph_func(params, t), T);

% Initial guess for parameters [D, B, A, b, C1, C2, omega_res1, omega_res2]
initial_guess = [1e-43, 1e-6, 1e-31, 1, 1e9, 1e10, 1e12, 4e12];

% Set bounds
lb = [1e-50, 1e-12, 1e-35, 1, 1e7, 1e8, 1e10, 3e12];
ub = [1e-40, 1e-3, 1e-28, 1000, 1e11, 1e12, 1e14, 5e12];

% Create opts structure
opts = optimoptions('lsqcurvefit', 'MaxFunctionEvaluations', 1e4, 'MaxIterations', 1e3);

% Use lsqcurvefit for the fitting
[fitted_params, resnorm, residual, exitflag, output] = lsqcurvefit(fun, initial_guess, T_data, k_ph_data, lb, ub, opts);

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Error using lsqncommon (line 35)
Objective function is returning Inf or NaN values at initial point. lsqcurvefit cannot continue.

Error in lsqcurvefit (line 322)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,optimgetFlag,caller,...

% Generate fitted curve (using log spacing)
T_fit = logspace(log10(min(T_data)), log10(max(T_data)), 100);
k_ph_fit = fun(fitted_params, T_fit);

[Aquatris]

You are dealing with really interesting magnitudes ranging from 1e-50 to 5e12. Checkout floating point precision to have a better understanding why that is not a good problem formulation to solve numerically, Unless you can change those I dont think you can get rid of the warnings. The default tolerance values generally lie around 1e-4 or so and since your values can be around 1e-50, you get warnings saying algorithm is reducing stuff.

Also in the end your code seems to find the parameters without errors, although I doubt it is an acceptable solution to the problem.

[Malay Agarwal]

Hi @Johnathan,

Could you please share the file "0T.txt" so that I can reproduce the issue?

[Johnathan]

Hi Malay,

Thank you for the consideration. Here is the file.

[Johnathan]

Hi @Malay Agarwal,

Sorry for not properly tagging you earlier, the file is above.

[Torsten]

@Johnathan

Are you sure T is in Kelvin and not in degreeC in your data file ? Experimental conditions of 2 K are hard to achieve :-)

[Johnathan]

Hi @Torsten

I am sure it is Kelvin, the data is from a cryogen we have in our lab from Quantum Design, we can get down to .35 K with LHe.

[Walter Roberson]

I experimented changing the integral() call to double(vpaintegral()) . That change avoided the warnings about nan or inf. It did, however, give a warning about problems with the finite difference in dimension 1, and another similar warning about the finite difference in dimension 3. The overall result in k_ph_fit was quite different than k_ph_data, showing only a simple exponential type curve, with no lobes at all.

[Johnathan]

Hi @Walter Roberson

Thank you for experimenting. I am unsure as to how this worked in python pretty well in comparison.

I think I got something similar to you by doing the same thing.

I am either getting the fit to come out to nearly 0 or around the 1e21 magnitude, with no dips in conductivity.

Here is the python code in case you may know of a better way of implementing it into MATLAB.

The below code should return this attached figure.

import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import pandas as pd

# Constants
k_B = 1.380649e-23  # Boltzmann constant (J/K)
hbar = 1.0545718e-34  # Reduced Planck constant (J·s)

# Parameters for FeCl2
theta_D = 235  # Debye temperature (K) fecl2, constant 
v_s = 4.817e3  # Sound velocity (m/s) fecl2, solved for 

# Parameters
L = 3.38e-4  # (m)
A = 7.03e-31  # (K^-1 s^2)
b = 3.2  # unitless
C1 = 3.029e+09  # s^-1 adjusts height/slope of second peak 
C2 = 1.548e+10  # s^-1 adjusts height/slope of overall conductivity 
omega_res1_0T = 3.5043973583e+12  # res1 for 0T fecl2 paper, constant
omega_res2_0T = 2.9114816436e+12  # res2 for 0T fecl2 paper, constant
D = 0.8e-43     # adjusts the higher-temperature peak
B = -5.298e-06  # adjusts the lowee_temperature peak


# Given omega values for 0T (in Hz)
omega_res1_0T = 3.5043973583e+12  # res1 for 0T fecl2 paper 
omega_res2_0T = 2.9114816436e+12  # res2 for 0T fecl2 paper

# Modify tau_tot_inv to use the given omega values for 0T
def tau_tot_inv(omega, T, H):
    tau_boundary_inv = v_s / L
    tau_defect_inv = D * omega**4
    tau_dislocation_inv = B * omega
    tau_umklapp_inv = A * T * omega**3 * np.exp(-theta_D / (b * T))
    
    tau_mag1_inv = C1 * (omega**4 / (omega**2 - omega_res1_0T**2)**2) * \
                   (np.exp(-hbar * omega_res1_0T / (k_B * T)) / 
                    (1 + np.exp(-hbar * omega_res1_0T / (k_B * T))))
    
    tau_mag2_inv = C2 * (omega**4 / (omega**2 - omega_res2_0T**2)**2) * \
                   (np.exp(-hbar * omega_res2_0T / (k_B * T)) / 
                    (1 + np.exp(-hbar * omega_res2_0T / (k_B * T))))
    
    return tau_boundary_inv + tau_defect_inv + tau_dislocation_inv + \
           tau_umklapp_inv + tau_mag1_inv + tau_mag2_inv

def integrand(x, T, H):
    omega = (x * k_B * T) / hbar
    F = (x**4 * np.exp(x)) / (np.exp(x) - 1)**2
    return F * (1 / tau_tot_inv(omega, T, H))

def k_ph(T, H):
    prefactor = (k_B / (2 * np.pi**2 * v_s)) * ((k_B * T / hbar)**3)
    integral, _ = integrate.quad(integrand, 0, theta_D / T, args=(T, H))
    return prefactor * integral

# Temperature range for calculation
temperatures = np.logspace(np.log10(1), np.log10(100), num=100)

# Calculate k_ph for each temperature (assuming H = 0 for now)
H = 0  # Magnetic field (T)
k_ph_values = [k_ph(T, H) for T in temperatures]




# Load data from the 0T.txt file, skipping the first row
data = np.loadtxt('0T.txt', skiprows=1)

# Extract temperature and thermal conductivity
temperature = data[:, 0]  # First column: Temperature (K)
thermal_conductivity = data[:, 1]  # Second column: \kappa_{xx} (W/m K)


# Plot both calculated and experimental data
plt.figure(figsize=(10, 6))
plt.plot(temperatures, k_ph_values, marker='', linestyle='-', color='b', label='Calculated')
plt.plot(temperature, thermal_conductivity, marker='o', linestyle='-', color='r', label='Raw')
plt.xscale('log')
# plt.yscale('log') 
plt.xlabel('T(K)')
plt.ylabel('$\kappa_{xx}$ (W/K m)')
plt.title('Data: Calculated vs Experimental')
plt.legend()
plt.grid(True)
plt.show()

[Torsten]

And what do you get in MATLAB if you set the parameters obtained from python as initial guesses ?

[Johnathan]

@Torsten

That is what I did originally, so the results were the same as what I got from the code in the original post. There definitely seems to be something wrong with the method itself as opposed to the guesses and bounds I believe.

[Johnathan]

Hi @Aquatris

That makes sense about the magnitudes and precision. Do you recommend a symbolic integrator instead?

Also, yes the values do print out but they are not reasonable for the fit, as it is quite off.

[Torsten]

@Johnathan

Please give us the solution for the parameters from your Python code for comparison that produced the figures you included.

[Johnathan]

Hi @Torsten

I just saw this comment now, sorry about that. Here is the python code for my closest fitting.

The figure created is attached as well.

import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import pandas as pd

# Constants
k_B = 1.380649e-23  # Boltzmann constant (J/K)
hbar = 1.0545718e-34  # Reduced Planck constant (J·s)

# Parameters for FeCl2
theta_D = 235  # Debye temperature (K) fecl2, constant 
v_s = 4.817e3  # Sound velocity (m/s) fecl2, solved for 

# Parameters
L = 3.38e-4  # (m)
A = 7.03e-31  # (K^-1 s^2)
b = 3.2  # unitless
C1 = 3.029e+09  # s^-1 adjusts height/slope of second peak 
C2 = 1.548e+10  # s^-1 adjusts height/slope of overall conductivity 
omega_res1_0T = 3.5043973583e+12  # res1 for 0T fecl2 paper, constant
omega_res2_0T = 2.9114816436e+12  # res2 for 0T fecl2 paper, constant
D = 0.8e-43     # adjusts the higher-temperature peak
B = -5.298e-06  # adjusts the lowee_temperature peak


# Given omega values for 0T (in Hz)
omega_res1_0T = 3.5043973583e+12  # res1 for 0T fecl2 paper 
omega_res2_0T = 2.9114816436e+12  # res2 for 0T fecl2 paper

# Modify tau_tot_inv to use the given omega values for 0T
def tau_tot_inv(omega, T, H):
    tau_boundary_inv = v_s / L
    tau_defect_inv = D * omega**4
    tau_dislocation_inv = B * omega
    tau_umklapp_inv = A * T * omega**3 * np.exp(-theta_D / (b * T))
    
    tau_mag1_inv = C1 * (omega**4 / (omega**2 - omega_res1_0T**2)**2) * \
                   (np.exp(-hbar * omega_res1_0T / (k_B * T)) / 
                    (1 + np.exp(-hbar * omega_res1_0T / (k_B * T))))
    
    tau_mag2_inv = C2 * (omega**4 / (omega**2 - omega_res2_0T**2)**2) * \
                   (np.exp(-hbar * omega_res2_0T / (k_B * T)) / 
                    (1 + np.exp(-hbar * omega_res2_0T / (k_B * T))))
    
    return tau_boundary_inv + tau_defect_inv + tau_dislocation_inv + \
           tau_umklapp_inv + tau_mag1_inv + tau_mag2_inv

def integrand(x, T, H):
    omega = (x * k_B * T) / hbar
    F = (x**4 * np.exp(x)) / (np.exp(x) - 1)**2
    return F * (1 / tau_tot_inv(omega, T, H))

def k_ph(T, H):
    prefactor = (k_B / (2 * np.pi**2 * v_s)) * ((k_B * T / hbar)**3)
    integral, _ = integrate.quad(integrand, 0, theta_D / T, args=(T, H))
    return prefactor * integral

# Temperature range for calculation
temperatures = np.logspace(np.log10(1), np.log10(100), num=100)

# Calculate k_ph for each temperature (assuming H = 0 for now)
H = 0  # Magnetic field (T)
k_ph_values = [k_ph(T, H) for T in temperatures]




# Load data from the 0T.txt file, skipping the first row
data = np.loadtxt('0T.txt', skiprows=1)

# Extract temperature and thermal conductivity
temperature = data[:, 0]  # First column: Temperature (K)
thermal_conductivity = data[:, 1]  # Second column: \kappa_{xx} (W/m K)


# Plot both calculated and experimental data
plt.figure(figsize=(10, 6))
plt.plot(temperatures, k_ph_values, marker='', linestyle='-', color='b', label='Calculated')
plt.plot(temperature, thermal_conductivity, marker='o', linestyle='-', color='r', label='Raw')
plt.xscale('log')
# plt.yscale('log') 
plt.xlabel('T(K)')
plt.ylabel('$\kappa_{xx}$ (W/K m)')
plt.title('Data: Calculated vs Experimental')
plt.legend()
plt.grid(True)
plt.show()

Answer 1

[Matt J]

I recommend using fminspleas, downloadable from,

https://www.mathworks.com/matlabcentral/fileexchange/10093-fminspleas

for this problem instead of lsqcurvefit.

This is more suitable for your problem in several ways. First, the model function depends nonlinearly on only 3 of your unknowns b, and , which fminspleas can exploit (once you make the change of variables C_3=1/L). Also, your function might not be continuous and differentiable in and , violating the assumptions of lsqcurvefit. In particular, your integrand has poles at and and it is not clear that the interval of integration avoids these poles.

Additionally, it will probably help with the numerical behavior of the exponential expressions to implement the objective function in terms of the transformed variables , and rather than in terms of b, and directly.

Obviously, you will then also need to transform the lb, ub bounds accordingly. I doubt you would ever want to allow any of the unknowns to go outside of [-15,+15]. The exp() outputs reach crazy values there.

[Johnathan]

@Matt J

Thank you for the suggestion.

My implementation is below

% Load data from 0T.txt
data = readmatrix('0T.txt', 'HeaderLines', 1);
T_data = data(:, 1); % First column: Temperature (K)
k_ph_data = data(:, 2); % Second column: \kappa_{xx} (W/mK)

% Constants
k_B = 1.380649e-23; % Boltzmann constant in J/K
hbar = 1.0545718e-34; % Reduced Planck constant in J·s
v_s = 4.817e3; % Sound velocity in m/s
theta_D = 235; % Debye temperature in K

% Define the integrand function with transformed variables
integrand = @(x, T, D, B, A, x0, C1, C2, C3, x1, x2) ...
    (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* ...
    (v_s * C3 + D*(k_B*T/hbar).^4.*x.^4 + B*(k_B*T/hbar).*x + ...
    A*T*(k_B*T/hbar).^3.*x.^3 .* exp(-x0./T) + ...
    C1 * (k_B*T/hbar).^4.*x.^4 ./ ((x.^2 - (x1*T/hbar).^2).^2 + 1e-10) .* ...
    exp(-x1./T) ./ (1 + exp(-x1./T)) + ...
    C2 * (k_B*T/hbar).^4.*x.^4 ./ ((x.^2 - (x2*T/hbar).^2).^2 + 1e-10) .* ...
    exp(-x2./T) ./ (1 + exp(-x2./T)));

% Define the k_ph function using integral
k_ph_func = @(params, T) (k_B / (2 * pi^2 * v_s)) * (k_B * T / hbar).^3 .* ...
    arrayfun(@(t) integral(@(x) integrand(x, t, params(1), params(2), params(3), ...
    params(4), params(5), params(6), params(7), params(8), params(9)), ...
    0, theta_D/t), T);

% Define the funlist for fminspleas
funlist = {@(params, T) k_ph_func(params, T)};

% Initial guess for parameters [D, B, A, x0, C1, C2, C3, x1, x2]
NLPstart = [0.8e-43, -5.298e-06, 7.03e-31, 73.4375, 3.029e+09, 1.548e+10, 2.9585798817e+03, 5.0128134748e-21, 4.1615451692e-21];

% Set bounds for the parameters
INLB = [-15, -15, -15, -15, -15, -15, -15, -15, -15];
INUB = [15, 15, 15, 15, 15, 15, 15, 15, 15];

% Define options for the optimizer
options = optimset('Display', 'iter');

% Run fitting using fminspleas
[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options);

% Combine INLP and ILP into fitted_params
fitted_params = [INLP; ILP];

% Generate fitted curve (using logarithmic spacing)
T_fit = logspace(log10(min(T_data)), log10(max(T_data)), 100);
k_ph_fit = k_ph_func(fitted_params, T_fit);

% Plot results
figure;
loglog(T_data, k_ph_data, 'o', 'DisplayName', 'Experimental Data');
hold on;
loglog(T_fit, k_ph_fit, '-', 'DisplayName', 'Fitted Curve');
xlabel('Temperature (K)');
ylabel('Thermal Conductivity \kappa_{ph} (W/m K)');
legend;
title('Temperature-Dependent Thermal Conductivity Fitting');

% Display fitted parameters
disp('Fitted Parameters:');
disp(fitted_params);

% Transform fitted parameters back to original variables
b = theta_D / fitted_params(4);
L = 1 / fitted_params(7);
omega_res1 = fitted_params(8) * k_B / hbar;
omega_res2 = fitted_params(9) * k_B / hbar;

disp('Transformed Parameters:');
disp(['b = ', num2str(b)]);
disp(['L = ', num2str(L)]);
disp(['omega_res1 = ', num2str(omega_res1)]);
disp(['omega_res2 = ', num2str(omega_res2)]);

The output and errors are

Unrecognized function or variable 'T'.

Error in untitled2>@(p)objective_function(p,T,k_ph) (line 39)

[params, resnorm] = fminsearch(@(p) objective_function(p, T, k_ph), initial_guess, options);

Error in fminsearch (line 209)

fv(:,1) = funfcn(x,varargin{:});

Error in untitled2 (line 39)

[params, resnorm] = fminsearch(@(p) objective_function(p, T, k_ph), initial_guess, options);

Unrecognized function or variable 'k_ph'.

Error in untitled2>@(p)objective_function(p,T_data,k_ph) (line 39)

[params, resnorm] = fminsearch(@(p) objective_function(p, T_data, k_ph), initial_guess, options);

Error in fminsearch (line 209)

fv(:,1) = funfcn(x,varargin{:});

Error in untitled2 (line 39)

[params, resnorm] = fminsearch(@(p) objective_function(p, T_data, k_ph), initial_guess, options);

Optimization terminated:

the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04

and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04

Error using vertcat

Dimensions of arrays being concatenated are not consistent.

Error in untitled2 (line 45)

fitted_params = [INLP; ILP];}

>>

[Matt J]

I don't understand how you've set up your funlist. You need to break up your integral into 6 terms so that your model looks like,

k_ph= A*f_1(INLP, T) + B*f_2(INLP, T) + D*f_3(INLP, T) + ...
      C_1*f_4(INLP, T) + C_2*f_5(INLP, T) + C_3*f_6(INLP, T)

where C_3=1/L and INLP is a vector of the remaining 3 parameters (transfomed as I recommended).

Once you have it in this form, the funlist shoudl be,

funlist={f_1, f_2,...,f_6}

Also, the INLB and INUB vectors should be of length 3 and specify bounds on the three INLP parameters only, not all nine. Similarly with NLPstart.

Then, after you've run fminspleas,

[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options);

the output INLP will contain the optimal solution for the three nonlinear paraters, while ILP will contain the 6 linear parameters [A,B,D,C_1,C_2,C_3].

[Sam Chak]

Hi @Johnathan

Would using a sufficiently high-order Taylor series to approximate the definite integral to a desired accuracy make the curve fitting relatively easier?

[Matt J]

@Sam Chak Probably not. The only terms of the integrand that are not already polynomials are rational functions, and rational functions are poorly approximated by polynomials, at least in the neighborhood of the poles.

[Johnathan]

Hi @Matt J

Sorry for the delay, here is my implementation based on your edits:

% Load data from 0T.txt
data = readmatrix('0T.txt', 'HeaderLines', 1);
T_data = data(:, 1); % First column: Temperature (K)
k_ph_data = data(:, 2); % Second column: \kappa_{xx} (W/mK)

% Constants
k_B = 1.380649e-23; % Boltzmann constant in J/K
hbar = 1.0545718e-34; % Reduced Planck constant in J·s
v_s = 4.817e3; % Sound velocity in m/s
theta_D = 235; % Debye temperature in K

% Define the 6 separate terms of the integrand
f1 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).^4.*x.^4;
f2 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).*x;
f3 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T*(k_B*T/hbar).^3.*x.^3 .* exp(-INLP(1)./T);
f4 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).^4.*x.^4 ./ ((x.^2 - (INLP(2)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(2)./T) ./ (1 + exp(-INLP(2)./T));
f5 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).^4.*x.^4 ./ ((x.^2 - (INLP(3)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(3)./T) ./ (1 + exp(-INLP(3)./T));
f6 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* v_s;

% Define the k_ph function using integral
k_ph_func = @(INLP, ILP, T) (k_B / (2 * pi^2 * v_s)) * (k_B * T / hbar).^3 .* ...
    arrayfun(@(t) ...
    ILP(1) * integral(@(x) f1(x, t, INLP), 0, theta_D/t) + ...
    ILP(2) * integral(@(x) f2(x, t, INLP), 0, theta_D/t) + ...
    ILP(3) * integral(@(x) f3(x, t, INLP), 0, theta_D/t) + ...
    ILP(4) * integral(@(x) f4(x, t, INLP), 0, theta_D/t) + ...
    ILP(5) * integral(@(x) f5(x, t, INLP), 0, theta_D/t) + ...
    ILP(6) * integral(@(x) f6(x, t, INLP), 0, theta_D/t), T);

% Define the funlist for fminspleas
funlist = {@(INLP, T) integral(@(x) f1(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f2(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f3(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f4(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f5(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f6(x, T, INLP), 0, theta_D/T)};

% Initial guess for nonlinear parameters [x0, x1, x2] (generated guesses)
NLPstart = [73.4375, 5.0128134748e-21, 4.1615451692e-21];

% Set bounds for the nonlinear parameters
INLB = [0, 0, 0];
INUB = [100, 1e-20, 1e-20];

% Define opts for the optimizer
options = optimset('Display', 'iter');

% fitting using fminspleas
[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options);

% Generate fitted curve (using log spacing)
T_fit = logspace(log10(min(T_data)), log10(max(T_data)), 100);
k_ph_fit = k_ph_func(INLP, ILP, T_fit);

% Plots
figure;
loglog(T_data, k_ph_data, 'o', 'DisplayName', 'Experimental Data');
hold on;
loglog(T_fit, k_ph_fit, '-', 'DisplayName', 'Fitted Curve');
xlabel('Temperature (K)');
ylabel('Thermal Conductivity \kappa_{ph} (W/m K)');
legend;
title('Temperature-Dependent Thermal Conductivity Fitting');

% Display fitted parameters
disp('Fitted Nonlinear Parameters (INLP):');
disp(INLP);
disp('Fitted Linear Parameters (ILP):');
disp(ILP);

% Transform fitted parameters back to original variables
x0 = INLP(1);
x1 = INLP(2);
x2 = INLP(3);
b = theta_D / x0;
L = 1 / ILP(6);
omega_res1 = x1 * k_B / hbar;
omega_res2 = x2 * k_B / hbar;

disp('Transformed Parameters:');
disp(['b = ', num2str(b)]);
disp(['L = ', num2str(L)]);
disp(['omega_res1 = ', num2str(omega_res1)]);
disp(['omega_res2 = ', num2str(omega_res2)]);

Which returns this:

Error using integral (line 84)

Limits of integration must be double or single scalars.

Error in untitled4>@(INLP,T)integral(@(x)f1(x,T,INLP),0,theta_D/T) (line 31)

funlist = {@(INLP, T) integral(@(x) f1(x, T, INLP), 0, theta_D/T), ...

Error in fminspleas/fminspleas_obj (line 264)

term = fi(INLP,xdata);

Error in fminsearch (line 209)

fv(:,1) = funfcn(x,varargin{:});

Error in fminspleas (line 293)

tINLP = fminsearch(@fminspleas_obj,tNLPstart,options);

Error in untitled4 (line 49)

[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options);

I tried a different way to alleviate some errors:

% Load data from 0T.txt
data = readmatrix('0T.txt', 'HeaderLines', 1);
T_data = data(:, 1); % First column: Temperature (K)
k_ph_data = data(:, 2); % Second column: \kappa_{xx} (W/mK)

% Constants
k_B = 1.380649e-23; % Boltzmann constant in J/K
hbar = 1.0545718e-34; % Reduced Planck constant in J·s
v_s = 4.817e3; % Sound velocity in m/s
theta_D = 235; % Debye temperature in K

% Define the 6 separate terms of the integrand
f1 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2);
f2 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* x.^(-3);
f3 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* x.^(-1) .* exp(-exp(INLP(1))./T);
f4 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) ./ ((x.^2 - (exp(INLP(2))*T/hbar).^2).^2 + 1e-10) .* exp(-exp(INLP(2))./T) ./ (1 + exp(-exp(INLP(2))./T));
f5 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) ./ ((x.^2 - (exp(INLP(3))*T/hbar).^2).^2 + 1e-10) .* exp(-exp(INLP(3))./T) ./ (1 + exp(-exp(INLP(3))./T));
f6 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2);

% funlist for fminspleas
funlist = {@(INLP, T) f1(theta_D/T, T, INLP), ...
           @(INLP, T) f2(theta_D/T, T, INLP), ...
           @(INLP, T) f3(theta_D/T, T, INLP), ...
           @(INLP, T) f4(theta_D/T, T, INLP), ...
           @(INLP, T) f5(theta_D/T, T, INLP), ...
           @(INLP, T) f6(theta_D/T, T, INLP)};

% k_ph function
k_ph_func = @(INLP, ILP, T) (k_B / (2 * pi^2 * v_s)) * (k_B * T / hbar).^3 .* ...
    (ILP(1) * funlist{1}(INLP, T) + ...
     ILP(2) * funlist{2}(INLP, T) + ...
     ILP(3) * funlist{3}(INLP, T) + ...
     ILP(4) * funlist{4}(INLP, T) + ...
     ILP(5) * funlist{5}(INLP, T) + ...
     (1/ILP(6)) * funlist{6}(INLP, T));

% Initial guess for nonlinear parameters [log(x0), log(x1), log(x2)]
% (generated guesses)
NLPstart = log([73.4375, 5.0128134748e-21, 4.1615451692e-21]);

% Set bounds for the nonlinear parameters
INLB = log([1e-6, 1e-30, 1e-30]);
INUB = log([1e3, 1e-10, 1e-10]);

% Define options for the optimizer
options = optimset('Display', 'iter');

% Wrapper function for fminspleas
wrapper = @(INLP, T) k_ph_func(INLP, [ones(1,5), 1], T);

% fitting using fminspleas
[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options, wrapper);

% Generate fitted curve (using log spacing)
T_fit = logspace(log10(min(T_data)), log10(max(T_data)), 100);
k_ph_fit = k_ph_func(INLP, ILP, T_fit);

% Plots
figure;
loglog(T_data, k_ph_data, 'o', 'DisplayName', 'Experimental Data');
hold on;
loglog(T_fit, k_ph_fit, '-', 'DisplayName', 'Fitted Curve');
xlabel('Temperature (K)');
ylabel('Thermal Conductivity \kappa_{ph} (W/m K)');
legend;
title('Temperature-Dependent Thermal Conductivity Fitting');

% Display fitted parameters
disp('Fitted Nonlinear Parameters (INLP):');
disp(exp(INLP));
disp('Fitted Linear Parameters (ILP):');
disp(ILP);

% Transform fitted parameters back to original variables
x0 = theta_D / exp(INLP(1));
x1 = exp(INLP(2)) * k_B / hbar;
x2 = exp(INLP(3)) * k_B / hbar;
b = exp(INLP(1));
L = 1 / ILP(6);
omega_res1 = x1;
omega_res2 = x2;

disp('Transformed Parameters:');
disp(['b = ', num2str(b)]);
disp(['L = ', num2str(L)]);
disp(['omega_res1 = ', num2str(omega_res1)]);
disp(['omega_res2 = ', num2str(omega_res2)]);

But i cannot suppress the errors:

Error using fminspleas

Too many input arguments.

Error in untitled5 (line 51)

[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options, wrapper);

I think I am misunderstanding something minor, but cannot quite figure it out. Thanks again for your help!

[Johnathan]

Hi @Sam Chak

As Matt says, I dont think the expansion would work for that portion of the function. But, if you can find a way to get a sufficiently high accuracy using such method, that would be of great help!

[Sam Chak]

Hi @Johnathan

Thank you for your update. The suggestion was purely theoretical, and I appreciate that @Matt J pointed out the issue of poor approximation of the rational function using Taylor series.

[Matt J]

@Johnathan

You seem to have invented a syntax that is not in the fminspleas documentation. In this line,

% fitting using fminspleas
[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options, wrapper);

you have given 9 input arguments, whereas minspleas only takes eight.

[Johnathan]

@Matt J

That was something I tried to alleviate the other errors, this what I tried before with the correct number of inputs.

% Load data from 0T.txt
data = readmatrix('0T.txt', 'HeaderLines', 1);
T_data = data(:, 1); % First column: Temperature (K)
k_ph_data = data(:, 2); % Second column: \kappa_{xx} (W/mK)

% Constants
k_B = 1.380649e-23; % Boltzmann constant in J/K
hbar = 1.0545718e-34; % Reduced Planck constant in J·s
v_s = 4.817e3; % Sound velocity in m/s
theta_D = 235; % Debye temperature in K

% Define the 6 separate terms of the integrand
f1 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).^4.*x.^4;
f2 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).*x;
f3 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T*(k_B*T/hbar).^3.*x.^3 .* exp(-INLP(1)./T);
f4 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).^4.*x.^4 ./ ((x.^2 - (INLP(2)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(2)./T) ./ (1 + exp(-INLP(2)./T));
f5 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).^4.*x.^4 ./ ((x.^2 - (INLP(3)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(3)./T) ./ (1 + exp(-INLP(3)./T));
f6 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* v_s;

% Define the k_ph function using integral
k_ph_func = @(INLP, ILP, T) (k_B / (2 * pi^2 * v_s)) * (k_B * T / hbar).^3 .* ...
    arrayfun(@(t) ...
    ILP(1) * integral(@(x) f1(x, t, INLP), 0, theta_D/t) + ...
    ILP(2) * integral(@(x) f2(x, t, INLP), 0, theta_D/t) + ...
    ILP(3) * integral(@(x) f3(x, t, INLP), 0, theta_D/t) + ...
    ILP(4) * integral(@(x) f4(x, t, INLP), 0, theta_D/t) + ...
    ILP(5) * integral(@(x) f5(x, t, INLP), 0, theta_D/t) + ...
    ILP(6) * integral(@(x) f6(x, t, INLP), 0, theta_D/t), T);

% Define the funlist for fminspleas
funlist = {@(INLP, T) integral(@(x) f1(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f2(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f3(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f4(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f5(x, T, INLP), 0, theta_D/T), ...
           @(INLP, T) integral(@(x) f6(x, T, INLP), 0, theta_D/T)};

% Initial guess for nonlinear parameters [x0, x1, x2]
NLPstart = [73.4375, 5.0128134748e-21, 4.1615451692e-21];

% Set bounds for the nonlinear parameters
INLB = [0, 0, 0];
INUB = [100, 1e-20, 1e-20];

% Define options for the optimizer
options = optimset('Display', 'iter');

% Run fitting using fminspleas
[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options);

% Generate fitted curve (using logarithmic spacing)
T_fit = logspace(log10(min(T_data)), log10(max(T_data)), 100);
k_ph_fit = k_ph_func(INLP, ILP, T_fit);

% Plot results
figure;
loglog(T_data, k_ph_data, 'o', 'DisplayName', 'Experimental Data');
hold on;
loglog(T_fit, k_ph_fit, '-', 'DisplayName', 'Fitted Curve');
xlabel('Temperature (K)');
ylabel('Thermal Conductivity \kappa_{ph} (W/m K)');
legend;
title('Temperature-Dependent Thermal Conductivity Fitting');

% Display fitted parameters
disp('Fitted Nonlinear Parameters (INLP):');
disp(INLP);
disp('Fitted Linear Parameters (ILP):');
disp(ILP);

% Transform fitted parameters back to original variables
x0 = INLP(1);
x1 = INLP(2);
x2 = INLP(3);
b = theta_D / x0;
L = 1 / ILP(6);
omega_res1 = x1 * k_B / hbar;
omega_res2 = x2 * k_B / hbar;

disp('Transformed Parameters:');
disp(['b = ', num2str(b)]);
disp(['L = ', num2str(L)]);
disp(['omega_res1 = ', num2str(omega_res1)]);
disp(['omega_res2 = ', num2str(omega_res2)]);

Which results in

Error using integral (line 84)

Limits of integration must be double or single scalars.

Error in untitled4>@(INLP,T)integral(@(x)f1(x,T,INLP),0,theta_D/T) (line 31)

funlist = {@(INLP, T) integral(@(x) f1(x, T, INLP), 0, theta_D/T), ...

Error in fminspleas/fminspleas_obj (line 264)

term = fi(INLP,xdata);

Error in fminsearch (line 209)

fv(:,1) = funfcn(x,varargin{:});

Error in fminspleas (line 293)

tINLP = fminsearch(@fminspleas_obj,tNLPstart,options);

Error in untitled4 (line 49)

[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options);

[Matt J]

Each funlist{i} needs to be modified with arrayfun the same way as you have done in kph_func. They all need to be able to accept a vector of T values.

It might be worth defining yourself a dedicated function for this that can be re-used

funlist = { vecIntegral(f1,theta_d), ...
            vecIntegral(f2,theta_d), ...
            vecIntegral(f3,theta_d), ...
            vecIntegral(f4,theta_d), ...
            vecIntegral(f5,theta_d), ...
            vecIntegral(f6,theta_d)};

function integratedf = vecIntegral(f,theta_d)

   integratedf = @(INLP,T) ...
          arrayfun(  @(t)integral(@(x) f(x, t, INLP), 0, theta_D/t) ,  T);

end

[Johnathan]

@Matt J

Her is my recent attempt:

% Load data from 0T.txt
data = readmatrix('0T.txt', 'HeaderLines', 1);
T_data = data(:, 1); % First column: Temperature (K)
k_ph_data = data(:, 2); % Second column: \kappa_{xx} (W/mK)

% Constants
k_B = 1.380649e-23; % Boltzmann constant in J/K
hbar = 1.0545718e-34; % Reduced Planck constant in J·s
v_s = 4.817e3; % Sound velocity in m/s
theta_D = 235; % Debye temperature in K

% Define the 6 separate terms of the integrand
f1 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).^4.*x.^4;
f2 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).*x;
f3 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T*(k_B*T/hbar).^3.*x.^3 .* exp(-INLP(1)./T);
f4 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).^4.*x.^4 ./ ((x.^2 - (INLP(2)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(2)./T) ./ (1 + exp(-INLP(2)./T));
f5 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* (k_B*T/hbar).^4.*x.^4 ./ ((x.^2 - (INLP(3)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(3)./T) ./ (1 + exp(-INLP(3)./T));
f6 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* v_s;

% Define the k_ph function using integral
k_ph_func = @(INLP, ILP, T) (k_B / (2 * pi^2 * v_s)) * (k_B * T / hbar).^3 .* ...
 arrayfun(@(t) ...
 ILP(1) * integral(@(x) f1(x, t, INLP), 0, theta_D/t) + ...
 ILP(2) * integral(@(x) f2(x, t, INLP), 0, theta_D/t) + ...
 ILP(3) * integral(@(x) f3(x, t, INLP), 0, theta_D/t) + ...
 ILP(4) * integral(@(x) f4(x, t, INLP), 0, theta_D/t) + ...
 ILP(5) * integral(@(x) f5(x, t, INLP), 0, theta_D/t) + ...
 ILP(6) * integral(@(x) f6(x, t, INLP), 0, theta_D/t), T);

% Define the funlist for fminspleas
funlist = {
    @(INLP, T) arrayfun(@(t) integral(@(x) f1(x, t, INLP), 0, theta_D/t), T), ...
    @(INLP, T) arrayfun(@(t) integral(@(x) f2(x, t, INLP), 0, theta_D/t), T), ...
    @(INLP, T) arrayfun(@(t) integral(@(x) f3(x, t, INLP), 0, theta_D/t), T), ...
    @(INLP, T) arrayfun(@(t) integral(@(x) f4(x, t, INLP), 0, theta_D/t), T), ...
    @(INLP, T) arrayfun(@(t) integral(@(x) f5(x, t, INLP), 0, theta_D/t), T), ...
    @(INLP, T) arrayfun(@(t) integral(@(x) f6(x, t, INLP), 0, theta_D/t), T)
};

% Initial guess for nonlinear parameters [x0, x1, x2]
NLPstart = [73.4375, 5.0128134748e-21, 4.1615451692e-21];

% Set bounds for the nonlinear parameters
INLB = [0, 0, 0];
INUB = [100, 1e-20, 1e-20];

% Define options for the optimizer
options = optimset('Display', 'iter');

% Run fitting using fminspleas
[INLP, ILP] = fminspleas(funlist, NLPstart, T_data, k_ph_data, INLB, INUB, [], options);

% Generate fitted curve (using logarithmic spacing)
T_fit = logspace(log10(min(T_data)), log10(max(T_data)), 100);
k_ph_fit = k_ph_func(INLP, ILP, T_fit);

% Plot results
figure;
loglog(T_data, k_ph_data, 'o', 'DisplayName', 'Experimental Data');
hold on;
loglog(T_fit, k_ph_fit, '-', 'DisplayName', 'Fitted Curve');
xlabel('Temperature (K)');
ylabel('Thermal Conductivity \kappa_{ph} (W/m K)');
legend;
title('Temperature-Dependent Thermal Conductivity Fitting');

% Display fitted parameters
disp('Fitted Nonlinear Parameters (INLP):');
disp(INLP);
disp('Fitted Linear Parameters (ILP):');
disp(ILP);

% Transform fitted parameters back to original variables
x0 = INLP(1);
x1 = INLP(2);
x2 = INLP(3);
b = theta_D / x0;
L = 1 / ILP(6);
omega_res1 = x1 * k_B / hbar;
omega_res2 = x2 * k_B / hbar;

disp('Transformed Parameters:');
disp(['b = ', num2str(b)]);
disp(['L = ', num2str(L)]);
disp(['omega_res1 = ', num2str(omega_res1)]);
disp(['omega_res2 = ', num2str(omega_res2)]);

and the result which produced values and a plot but not accurate

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Iteration Func-count f(x) Procedure

0 1 2064.85

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

1 4 2064.4 initial simplex

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

2 5 2064.4 reflect

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

3 7 2064.1 expand

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

4 9 2063.37 expand

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

5 10 2063.37 reflect

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

6 12 2062.1 expand

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

7 14 2060.79 expand

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

8 15 2060.79 reflect

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

9 17 2057.7 expand

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

10 19 2056.04 expand

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

11 21 2054.49 expand

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

12 23 2053.98 reflect

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

13 24 2053.98 reflect

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

14 26 2053.98 contract outside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

15 28 2053.98 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

16 30 2053.98 contract outside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

17 32 2053.98 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

18 33 2053.98 reflect

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

19 35 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

20 37 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

21 39 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

22 41 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

23 43 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

24 45 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

25 47 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

26 49 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

27 51 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

28 53 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

29 55 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

30 57 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

31 59 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

32 61 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

33 63 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

34 65 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

35 67 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

36 69 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

37 71 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

38 73 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

39 75 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

40 77 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

41 79 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

42 81 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

43 83 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

44 85 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

45 87 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

46 89 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

47 91 2053.97 contract inside

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

48 93 2053.97 contract inside

Optimization terminated:

the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04

and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04

Warning: Rank deficient, rank = 2, tol = 3.127902e+41.

Fitted Nonlinear Parameters (INLP):

100.0000 0.0000 0.0000

Fitted Linear Parameters (ILP):

1.0e-41 *

0.0000

0

0.1186

0

0

0

Transformed Parameters:

b = 2.35

L = Inf

omega_res1 = 6.5454e-10

omega_res2 = 6.8216e-10

[Matt J]

I wonder if your bounds and NLPstart take into account the change of variables you made.

INLB = [0, 0, 0];
INUB = [100, 1e-20, 1e-20];

With such small bounds on INLP(2) and INLP(3), the exponentials where they appear will basically be 1 all the time.

[Johnathan]

@Matt J

I get the same warinings and very similar fit if i set the bounds as:

% Set bounds for the nonlinear parameters
INLB = [0, 0, 0];
INUB = [100, 1, 1];

[Matt J]

I notice now that your integrand is undefined at x=0 (which is also the left limit of integration), so the integral is improper. You'll have to reformulate it somehow, perhaps by using a Taylor series approximation to the exponentials near x=0. I believe this leads to,

[Matt J]

Also, maybe get rid of the factors of (k_B/hbar).^n? All these are doing are rescaling the linear coefficient ILP with vastly different orders of magnitude.

f1 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.^4.*x.^4;
f2 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.*x;
f3 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.^4.*x.^3 .* exp(-INLP(1)./T);
f4 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.^4.*x.^4 ./ ((x.^2 - (INLP(2)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(2)./T) ./ (1 + exp(-INLP(2)./T));
f5 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.^4.*x.^4 ./ ((x.^2 - (INLP(3)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(3)./T) ./ (1 + exp(-INLP(3)./T));
f6 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* v_s;

[Matt J]

Also, you currently have a lower bound for omega_res1, omega_res2 of 0. Don't they tend to be strictly greater than theta_D so that you won't be integrating through the singularities at omega_res1, omega_res2?

[Torsten]

According to the original question, the integration variable x equals h*omega/(k_B*T). This would make it necessary to rewrite the complete integrand in x instead of T. But maybe I'm mistaken.

[Matt J]

Yes, you're right, so INLP(2:3) should be >Theta_D

[Johnathan]

@Torsten

Do you mean with respect to this? Having to rewrite it all in terms of x.

f1 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.^4.*x.^4;
f2 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.*x;
f3 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.^4.*x.^3 .* exp(-INLP(1)./T);
f4 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.^4.*x.^4 ./ ((x.^2 - (INLP(2)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(2)./T) ./ (1 + exp(-INLP(2)./T));
f5 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* T.^4.*x.^4 ./ ((x.^2 - (INLP(3)*T/hbar).^2).^2 + 1e-10) .* exp(-INLP(3)./T) ./ (1 + exp(-INLP(3)./T));
f6 = @(x, T, INLP) (x.^4 .* exp(x) ./ (exp(x) - 1).^2) .* v_s;

[Torsten]

Usually, all the T's in your integrand should be replaced by h*omega/(k_B*x), shouldn't they ? But the upper limit of the integrand (Theta_D/T) also depends on T (which varies with x). Strange ...

[Matt J]

No, x is the variable of integration and T is constant with respect to the integration. It is omega that needs to be replaced with k_B*T*x/h.