Cannot figure out how to fit a complicated custom equation.
120 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1758559/image.png)
where
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1758564/image.png)
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);
% 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);
12 commentaires
Torsten
le 24 Août 2024 à 13:38
Please give us the solution for the parameters from your Python code for comparison that produced the figures you included.
Réponses (1)
Matt J
le 24 Août 2024 à 0:13
Modifié(e) : Matt J
le 24 Août 2024 à 1:13
I recommend using fminspleas, downloadable from,
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759014/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759019/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759014/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759019/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759014/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759019/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759064/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759069/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759074/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759014/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759019/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1759079/image.png)
4 commentaires
Sam Chak
le 24 Août 2024 à 7:20
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?
Voir également
Catégories
En savoir plus sur Eigenvalue Problems dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!