Fitting an anonymous function: Debye Callaway Model
18 views (last 30 days)
Show older comments
Hello,
I'm trying to fit the Debye Callaway Model of thermal conductivity-and relatively new at Matlab. I have the temperature and thermal conductivity, and have to find the two fitting parameters A and B, and I keep getting the following errors. Can someone please help me out how to better build my code?
TC is my excel table which has the temperature and Thermal conductivity data. Col 1 is Temperature and Column 6 is Thermal conductivity.
Thank you so much.
clear;
close all;
clc;
%Constants and parameters for calculation
TC=readmatrix('ExperimentalGe_Thermal Modeling.xlsx');
TD = 360; % Debye temperature in Kelvin
soundspeed = 4983; % Speed of sound in m/s
grain = 10e-5; % Grain size in meters
Planck = 1.05457e-34;
BoltzConst = 1.38066e-23;
C= BoltzConst^4 / (2*(pi^2).*Planck.^3*soundspeed);%Constants in the equation clubbed together
omega= BoltzConst.*TC(:,1)./Planck;%phonon frequency
%Klemen's Scattering Model
Bound = soundspeed/grain;
%PD = @(x) A.* omega.^4 *x.^4;
%Ump = @(x) B.*omega.^2 * T * x.^2 * exp(-TD/(3.0*TC(:,1)));
%totalTau = @(x) (Bound + Ump(x)+ PD(x)).^(-1);
%Fittype
FT=fittype(@(A,B,x,T)C.*T.^3.*(integral(@(x)(((Bound + B.*omega.^2 * T * x.^2 * exp(-TD/(3.0*T))+ A.* omega.^4 *x.^4).^(-1).*x.^4.*exp(x))./(exp(x)- 1).^2),0,TD./T)),'coefficients',{'A','B'},'independent',{'x','T'});
%Curvefitting
f=fit(TC(:,1),TC(:,6),FT )
Errors:
Error using fittype/testCustomModelEvaluation
Expression
C.*T.^3.*(integral(@(x)(((Bound+B.*omega.^2*T*x.^2*exp(-TD/(3.0*T))+A.*omega.^4*x.^4).^(-1).*x.^4.*exp(x))./(exp(x)-1).^2),0,TD./T)) is
not a valid MATLAB expression, has non-scalar coefficients, or cannot be evaluated:
Error while trying to evaluate FITTYPE function :
Limits of integration must be double or single scalars.
Error in fittype>iCreateFittype (line 373)
testCustomModelEvaluation( obj );
Error in fittype (line 330)
obj = iCreateFittype( obj, varargin{:} );
Error in TCGe (line 21)
FT=fittype(@(A,B,x,T)C.*T.^3.*(integral(@(x)(((Bound + B.*omega.^2 * T * x.^2 * exp(-TD/(3.0*T))+ A.* omega.^4 *x.^4).^(-1).*x.^4.*exp(x))./(exp(x)- 1).^2),0,TD./T)),'coefficients',{'A','B'},'independent',{'x','T'});
Caused by:
Error using fittype/evaluate>I_ERROR_FCN_
Error while trying to evaluate FITTYPE function :
Limits of integration must be double or single scalars.
>>
0 Comments
Answers (1)
Star Strider
on 15 Apr 2024
I do not have your data so I cannot run your code.
However it is always appropriate to compltely vectorise objective function operations (array or element-wise operations) unless you specifically intend some operations to be matrix operations —
C.*T.^3.*(integral(@(x)(((Bound+B.*omega.^2.*T.*x.^2.*exp(-TD./(3.0*T))+A.*omega.^4*x.^4).^(-1).*x.^4.*exp(x))./(exp(x)-1).^2),0,TD./T))
Some were not vectorised.
Also, rather than raising an expression or sub-expression to the (-1) power, divide it instead:
C.*T.^3.*(integral(@(x)((1./(Bound+B.*omega.^2.*T.*x.^2.*exp(-TD./(3.0*T))+A.*omega.^4*x.^4).*x.^4.*exp(x))./(exp(x)-1).^2),0,TD./T))
Check that last expression to be certain it is correct. The ‘1./(expression)’ is more computationally efficient than ‘(expression).^(-1)’.
.
6 Comments
Star Strider
on 16 Apr 2024
The fundamental problem that I found is the second limit of integration, that being ‘TD./T’. Since ‘T’ is a vector, that is going to cause problems, since the limits must by definition be scalars. Using arrayfun to integrate it with respect to each element of ‘T’ works, however it returns a (63x63) matrix when it should only return a (63x1) vector in each iteration of the objective function. (I experimented with fitnlm because I don’t have sufficient familiarity with the Curve Fitting Toolbox functions and syntax to work with them efficiently. I understand fitnlm well enough.)
I am posting my results to describe the problem. There may be a way to make this work, however it is late here (UCT-6) so I am calling it a day.
clear;
close all;
clc;
%Constants and parameters for calculation
% TC=readmatrix('ExperimentalGe_Thermal Modeling.xlsx');
TC = readmatrix('TCvsT.xlsx')
TD = 360; % Debye temperature in Kelvin
soundspeed = 4983; % Speed of sound in m/s
grain = 10e-5; % Grain size in meters
Planck = 1.05457e-34;
BoltzConst = 1.38066e-23;
C= BoltzConst^4 / (2*(pi^2).*Planck.^3*soundspeed);%Constants in the equation clubbed together
omega= BoltzConst.*TC(:,1)./Planck;%phonon frequency
%Klemen's Scattering Model
Bound = soundspeed/grain;
%PD = @(x) A.* omega.^4 *x.^4;
%Ump = @(x) B.*omega.^2 * T * x.^2 * exp(-TD/(3.0*TC(:,1)));
%totalTau = @(x) (Bound + Ump(x)+ PD(x)).^(-1);
% objfcn = @(A,B,T)C.*T.^3.*cell2mat(arrayfun(@(t)integral(@(x)((1./(Bound + B.*omega.^2 .* t * x.^2 .* exp(-TD./(3.0*t))+ A.* omega.^4 *x.^4).*x.^4.*exp(x))./(exp(x)- 1).^2),0,TD./t, 'ArrayValued',1),T, 'Unif',0))
objfcn = @(A,B,T)cell2mat(arrayfun(@(t)C.*t.^3.*integral(@(x)((1./(Bound + B.*omega.^2 .* t * x.^2 .* exp(-TD./(3.0*t))+ A.* omega.^4 *x.^4).*x.^4.*exp(x))./(exp(x)- 1).^2),0,TD./t, 'ArrayValued',1),T, 'Unif',0).')
fcn = @(b,T) objfcn(b(1),b(2),T);
Q = fcn(rand(2,1),TC(:,2))
mdl = fitnlm(TC(:,1),TC(:,2),fcn,rand(2,1))
return
%Fittype
FT=fittype(@(A,B,T)C.*T.^3.*(integral(@(x)((1./(Bound + B.*omega.^2 .* T * x.^2 .* exp(-TD./(3.0*T))+ A.* omega.^4 *x.^4).*x.^4.*exp(x))./(exp(x)- 1).^2),0,TD./T)),'coefficients',{'A','B'},'independent',{'T'});
%Curvefitting
f=fit(TC(:,1),TC(:,2),FT )
.
See Also
Categories
Find more on Linear and Nonlinear Regression in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!