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.

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.