A lei de arrefecimento de Newton fracionária

A lei de arrefecimento de Newton estabelece que a variação de Temperatura ao longo do tempo de um material ao esfriar é proporcional a diferença de temperatura entre o material e a vizinhança, isto é 

$\frac{dT(t)}{dt}=\lambda (T(t)-T_a)$

A solução desta equação diferencial pode ser encontrada pelo método operacional usando transformada de Laplace, 

O mesmo método pode ser usado para encontrar a solução da lei de arrefecimento de Newton fracionária

$\frac{d^\alpha T(t)}{dt^\alpha}=\lambda^\alpha (T(t)-T_a)$

onde $0<\alpha<1$,

onde $E_{\alpha,1}$ é a função de Mittag-Leffler. Segue abaixo o programa que faz um ajuste dados experimentais obtidos com auxílio de uma placa Arduino$(c)$ pelas equação de ordem fracionária. Os dados experimentais foram obtidos por Higor Monteiro para o seu TCC.

Um  resultado inicial é apresentado na figura abaixo $($com $\lambda=1.1160\times 10^{-4}$ e $\alpha=0.9939$$)$. Resultados usando dados parciais e cálculo aproximado da função de Mittag-Leffler não mostraram melhorias significativas. Uma discussão mais aprofundada da questão pode ser encontrada na referência JOURNAL OF APPLIED PHYSICS 123, 064901 (2018). Em breve, irei publicar mais detalhes sobre o Cálculo Fracionário. 

function fracnewton
global t T T0 Ta

load -ascii HomeLRN.txt HomeLRN
tdado=HomeLRN(:,1);
mdado=HomeLRN(:,2);
Tdado=HomeLRN(:,3);
Tadado=HomeLRN(:,4);
Udado=HomeLRN(:,5);

inix=find(Tdado==max(Tdado));
ini=max(inix);

r=400;
t0=tdado(ini);
t=tdado(ini:r:end)-t0;
%length(t)
T=Tdado(ini:r:end);
T0=Tdado(ini);
Ta=mean(Tadado(ini:r:end));
m=mdado(ini:r:end);
U=Udado(ini:r:end);

par0=[1e-4,0.9]; %valor inicial dos parâmetros

LB=[1e-5,0];
UB=[1e-3,1];

options=optimoptions('fmincon','Algorithm','interior-point','TolFun', 1e-10, 'TolX', 1e-12);
[par,Feval] = fmincon(@(x) myfun(x),par0,[],[],[],[],LB,UB,[],options);

lambda=par(1)
alpha=par(2)
Feval

Tcala=(T0-Ta)*mlf(alpha,1,-(lambda^alpha)*t.^alpha,12)+Ta;

figure(1)
plot(t,Tcala,'k-',t,T,'ko')
xlabel('t')
ylabel('T')

end

function E=myfun(z)
global t T T0 Ta
lambda=z(1);
alpha=z(2);
Tcala=(T0-Ta)*mlf(alpha,1,-(lambda^alpha)*t.^alpha,12)+Ta;
E=sqrt(sum(((Tcala-T)./T).^2))
end