Fitting an anonymous function: Debye Callaway Model

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.

>>

Answer 1

[Star Strider]

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)’.

.

[PJ]

I corrected the code to this:

FT = fittype(@(A,B,x,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',{'x','T'});

However, I'm not getting the error-

Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.

When I add another Parenthesis before Co-effecients, I end up getting these errors:

Error using fittype/testCustomModelEvaluation

Expression

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)) 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) ((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',{'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.

[Torsten]

x is a formal integration variable. It cannot be chosen as the independent variable of your model equation because it is integrated out and does not appear in the result for the expression

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)

[Star Strider]

@Torsten — Thank you!

I should have picked up on that, and didn’t. (I was concentrating on the syntax.)

We still don’t have the data.

[PJ]

@Star StriderAt

[PJ]

@Star Strider@Torsten I ran both the corrections, there's some delimiter error, I will correct that on the big screen tomorrow. Unable to find where I'm missing the paranthesis.

[Star Strider]

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')

TC = 63x2

317.3372 12.9288 315.0257 12.9096 312.7226 12.9003 310.4205 12.8891 308.1222 12.8786 305.8313 12.8726 303.5439 12.8638 301.2533 12.8590 298.9666 12.8490 296.6703 12.8363

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

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).')

objfcn = function_handle with value:
    @(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))

Q = 63x63

1.0e-38 * 0.0215 0.0214 0.0213 0.0212 0.0212 0.0211 0.0211 0.0210 0.0210 0.0209 0.0208 0.0207 0.0207 0.0206 0.0205 0.0204 0.0203 0.0203 0.0202 0.0202 0.0202 0.0201 0.0200 0.0200 0.0199 0.0198 0.0197 0.0195 0.0194 0.0192 0.0221 0.0220 0.0219 0.0219 0.0218 0.0217 0.0217 0.0216 0.0216 0.0215 0.0214 0.0213 0.0213 0.0212 0.0211 0.0210 0.0209 0.0209 0.0208 0.0208 0.0208 0.0207 0.0206 0.0205 0.0205 0.0204 0.0202 0.0201 0.0199 0.0198 0.0228 0.0226 0.0226 0.0225 0.0224 0.0224 0.0223 0.0223 0.0222 0.0221 0.0221 0.0220 0.0219 0.0218 0.0218 0.0217 0.0216 0.0215 0.0215 0.0214 0.0214 0.0213 0.0213 0.0212 0.0211 0.0210 0.0208 0.0207 0.0205 0.0204 0.0235 0.0233 0.0233 0.0232 0.0231 0.0231 0.0230 0.0230 0.0229 0.0228 0.0227 0.0226 0.0226 0.0225 0.0224 0.0223 0.0222 0.0222 0.0221 0.0221 0.0220 0.0219 0.0219 0.0218 0.0217 0.0216 0.0215 0.0213 0.0211 0.0210 0.0242 0.0240 0.0240 0.0239 0.0238 0.0238 0.0237 0.0237 0.0236 0.0235 0.0234 0.0233 0.0232 0.0232 0.0231 0.0230 0.0229 0.0228 0.0228 0.0227 0.0227 0.0226 0.0225 0.0225 0.0224 0.0223 0.0221 0.0220 0.0218 0.0216 0.0249 0.0248 0.0247 0.0246 0.0245 0.0245 0.0244 0.0244 0.0243 0.0242 0.0241 0.0240 0.0239 0.0239 0.0238 0.0237 0.0236 0.0235 0.0235 0.0234 0.0234 0.0233 0.0232 0.0231 0.0231 0.0229 0.0228 0.0226 0.0224 0.0222 0.0257 0.0255 0.0254 0.0253 0.0253 0.0252 0.0252 0.0251 0.0250 0.0249 0.0249 0.0248 0.0247 0.0246 0.0245 0.0244 0.0243 0.0242 0.0242 0.0241 0.0241 0.0240 0.0239 0.0238 0.0238 0.0236 0.0235 0.0233 0.0231 0.0229 0.0265 0.0263 0.0262 0.0261 0.0260 0.0260 0.0259 0.0259 0.0258 0.0257 0.0256 0.0255 0.0254 0.0254 0.0253 0.0252 0.0250 0.0250 0.0249 0.0249 0.0248 0.0247 0.0247 0.0246 0.0245 0.0244 0.0242 0.0240 0.0238 0.0236 0.0273 0.0271 0.0270 0.0269 0.0269 0.0268 0.0267 0.0267 0.0266 0.0265 0.0264 0.0263 0.0262 0.0261 0.0261 0.0259 0.0258 0.0258 0.0257 0.0257 0.0256 0.0255 0.0254 0.0253 0.0253 0.0251 0.0250 0.0248 0.0246 0.0244 0.0281 0.0280 0.0279 0.0278 0.0277 0.0276 0.0276 0.0275 0.0274 0.0273 0.0272 0.0271 0.0270 0.0270 0.0269 0.0268 0.0266 0.0266 0.0265 0.0265 0.0264 0.0263 0.0262 0.0261 0.0260 0.0259 0.0257 0.0255 0.0253 0.0251

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

mdl = fitnlm(TC(:,1),TC(:,2),fcn,rand(2,1))

Error using nlinfit (line 219)
MODELFUN must be a function that returns a vector of fitted values the same size as Y (63-by-1). The model function you provided returned a result that was 63-by-63.

One common reason for a size mismatch is using matrix operators (*, /, ^) in your function instead of the corresponding elementwise operators (.*, ./, .^).

Error in NonLinearModel/fitter (line 1168)
nlinfit(X,y,F,b0,opts,wtargs{:},errormodelargs{:});

Error in classreg.regr.FitObject/doFit (line 94)
model = fitter(model);

Error in NonLinearModel.fit (line 1487)
model = doFit(model);

Error in fitnlm (line 99)
model = NonLinearModel.fit(X,varargin{:});

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 )

.