Fitting an anonymous function: Debye Callaway Model

조회 수: 11 (최근 30일)
PJ 2024년 4월 15일
댓글: Star Strider 2024년 4월 16일
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
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.
>>

댓글을 달려면 로그인하십시오.

답변 (1개)

Star Strider 2024년 4월 15일
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이전 댓글 4개 표시이전 댓글 4개 숨기기
PJ 2024년 4월 16일
@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 2024년 4월 16일
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 = 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 )
.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Linear and Nonlinear Regression에 대해 자세히 알아보기

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by