Curve fitting with complexe custom function

2 views (last 30 days)
I want to fit some experimental data to a custom equation of this type:
Where a(T) and c(T) are some complicated functions of T (T is my x-axis parameter).
I calculate a(T) from a for loop. In my code a(T) and c(T) are represented by preExpPu(i,j,k) and preExpU(i,j,k). So basically the equation from the Figure above is my equation represented in the for loop by condTot(i,j,k). I know the values of a(T) and c(T) because I have calculated them. However, what I would like to fit is the values of b and d.... (in the exponential functions). The difficult part is that a(T) and c(T) are functions of T so I could not fit the whole data to a function of a type (a/T) exp(-b/T) + (c/T) exp(-d/T) because a and b are simply not constants, by vary with T... Is there any way to fit directly in the for loop do get the values of b and d?
for i=1:numel(T)
for j=1:numel(x)
for k=1:numel(y)
K(i)=exp((-H/(k1*T(i)))+1);
%Calcul de la fréquence moyenne des phonons (s-1)
v(j,k)=((2*pi*k2)/h)*TempDeb(j,k);
%Calcul des concentrations en trous/electrons
p(i,j,k)=-x(j)+(K(i)-sqrt(4*K(i)*x(j)*(K(i)*(x(j)-2*y(k)+1)+2*(y(k)-x(j))-1)+4*K(i)*y(k)*((y(k)-1)*(K(i)-1))+(K(i)^2)+4*(x(j)^2)))/(2*(K(i)-1));
n(i,j,k)=p(i,j,k)+2*x(j);
% % %Calcul des conductivité électriques dues à Pu et à U (S/m)
preExpPu(i,j,k)=(((4*(e^2)*v(j,k))/(a(j,k)*k2))*((n(i,j,k))/(y(k)))*(1-((n(i,j,k))/(y(k)))));
preExpU(i,j,k)=(((4*(e^2)*v(j,k))/(a(j,k)*k2))*((p(i,j,k))/(1-y(k)))*(1-((p(i,j,k))/(1-y(k)))));
condPu(i,j,k)= (preExpPu(i,j,k))/(T(:,i))*exp(-c/(k1*(T(:,i))));
condU(i,j,k)= (preExpU(i,j,k))/(T(:,i))*exp(-d/(k1*(T(:,i))));
% % %Calcul de la conductivité électrique totale (S/m)
condTot(i,j,k)=condPu(i,j,k)+condU(i,j,k);
end
end
end
  2 Comments
Plamen Bonev
Plamen Bonev on 29 Jun 2021
x=[0, 0.01, 0.03, 0.04, 0.05, 0.07];
% x=[0.08, 0.09, 0.1, 0.15, 0.16, 0.17];
% x=[1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2];
T=[500 : 50 : 3223]
y=[0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9];
Na=6.02e23;
R=8.314;
%Enthalpie de réaction en eV
H=1.5;
%h en J.s
h=6.63e-34
%k en eV/K
k1=8.6173303e-5;
%k en J/K
k2=1.38e-23;
%Charge élémentaire en Coulombs
e=1.6e-19;
%Energies de migration en eV
EhPu=0.8;
EhU=0.3;
%rayons en Angstrom (m)
rU4=0.1001e-9;
rPu4=0.096e-9;
rPu3=0.11e-9;
rO2=0.1368e-9;
rOv=0.1367e-9;
%Masses atomiques en Amu
Mu=238;
Mpu=244;
Mo=15.9;
eps=25.85;
format shortE
p=zeros(length(T),length(x),length(y));
n=zeros(length(T),length(x),length(y));
K=zeros(length(T),1);
condPu=zeros(length(T),length(x),length(y));
condU=zeros(length(T),length(x),length(y));
condTot=zeros(length(T),length(x),length(y));
condElectron=zeros(length(T),length(x),length(y));
yU4=zeros(length(y),1);
yPu4=zeros(length(x), length(y));
yPu4stoech=zeros(length(y),1);
yPu3=zeros(length(x),1);
yO2=zeros(length(x),1);
yO2stoech=(2/3);
yOv=zeros(length(x),1);
rmeanCat=zeros(length(x), length(y));
rmeanAn=zeros(length(x),1);
a=zeros(length(x),length(y));
vl=zeros(length(x),length(y));
vt=zeros(length(x),length(y));
N=zeros(length(x),1);
TempDeb=zeros(length(x),length(y));
v=zeros(length(x),length(y));
VaMetre=zeros(length(x),length(y));
VaAngstr=zeros(length(x),length(y));
Mmoy=zeros(length(y),1);
rmoy=zeros(length(y),1);
SommeGamma=zeros(length(x),length(y));
V=zeros(length(x),length(y));
vp=zeros(length(x),length(y));
C=zeros(length(x),length(y));
A=zeros(length(x),length(y));
LTE=zeros(length(x),length(y));
a0=zeros(length(x),length(y));
a1=zeros(length(x),length(y));
a2=zeros(length(x),length(y));
a3=zeros(length(x),length(y));
alpha=zeros(length(x),length(y));
MmoyMOXkG=zeros(length(x),length(y));
MmoyMOXamu=zeros(length(x),length(y));
rhoTh=zeros(length(x),length(y));
BulkModulus=zeros(length(x),length(y));
Cv=zeros(length(x),1);
Vm=zeros(length(x),length(y));
Grun=zeros(length(x),length(y));
constB=zeros(length(x),length(y));
B=zeros(length(x),length(y));
CondPhonique=zeros(length(T),length(x),length(y));
preExpPu=zeros(length(T),length(x),length(y));
preExpU=zeros(length(T),length(x),length(y));
% %y=0.2
for i=1:numel(T)
for j=1:numel(x)
for k=1:numel(y)
%Pour le calcul du paramètre de maille
yU4(k)=(1-y(k));
yPu4(j,k)=(y(k)-2*x(j));
yPu4stoech(k)=y(k)/3;
yPu3(j)=(2*x(j));
yO2(j)=(2-x(j))/2;
yOv(j)=x(j)/2;
rmeanCat(j,k)=yPu3(j)*rPu3+yPu4(j,k)*rPu4+yU4(k)*rU4;
rmeanAn(j)=yO2(j)*rO2;
%Pour le calcul de Gamma(i)
rmoy(k)=(yU4(k)/3)*rU4+yPu4stoech(k)*rPu4+(yO2stoech)*rO2;
Mmoy(k)=(yU4(k)/3)*Mu+yPu4stoech(k)*Mpu+(yO2stoech)*Mo;
%Pour le calcul de la fréquence moyenne des phonons
vl(j,k)=5358*(1-0.7279*x(j))*(1+0.040*y(k));
vt(j,k)=2750*(1-1.0545*x(j))*(1+0.043*y(k));
N(j)=4*(3-x(j));
%Calcul du paramètre de maille
a(j,k)=(4/(3^(1/2)))*(rmeanCat(j,k)+rmeanAn(j));
% Calcul du volume atomique moyen
VaMetre(j,k)=((a(j,k))^3)/N(j); %en m^3
VaAngstr(j,k)=((a(j,k)*10^(10))^3)/N(j); %en A^3
%Calcul du volume atomique moyen par molecule de MOX
V(j,k)=((a(j,k))^3)/4;
% Calcul de la somme de Gamma(i);
SommeGamma(j,k)=(((x(j))/3)*(((-Mo^2)/(Mmoy(k))^2)+(eps/(rmoy(k)^2))*(2*(rPu3^2-rPu4^2)-(rOv^2-rO2^2))))+(0.035/0.886);
%
%Calcul de la temperature de Debye (K)
TempDeb(j,k)=(h/k2)*((9*N(j))/(4*pi*(a(j,k)^3)*(1/(vl(j,k))^3+2/(vt(j,k))^3)))^(1/3);
%Calcul de la vitesse moyenne des phonons (m/s)
vp(j,k)=((2*pi*k2*TempDeb(j,k))/(h))*((V(j,k))/(6*(pi^2)))^(1/3);
%Calcul de C
C(j,k)=((pi^2)*VaMetre(j,k)*TempDeb(j,k))/(3*(vp(j,k)^2)*h);
%Calcul de A
A(j,k)=C(j,k)*SommeGamma(j,k);
%--------------------------------------------------------------------------
%Pour le calcul du LTE
a0(j,k)=(-2.8809e-3)+((0.0301e-3)*y(k))-(4.3954e-3)*x(j)+(0.0156e-3*((y(k))^2))-(15.1759e-3)*((x(j))^2)+(2.5642e-3)*y(k)*x(j);
a1(j,k)=(9.5024e-6)+((-0.1864e-6)*y(k))+((15.8173e-6)*x(j))+((-0.0229e-6)*(y(k))^2)+((7.6258e-6)*(x(j))^2)+((-7.5789e-6)*y(k)*x(j));
a2(j,k)=(2.0894e-10)+(2.9483e-10)*y(k)+(-19.9227e-10)*x(j)+(-1.0355e-10)*(y(k))^2+(73.8931e-10)*(x(j))^2+(11.6442e-10)*y(k)*x(j);
a3(j,k)=(4.4096e-13)+(-1.4263e-13)*y(k)+(23.5638e-13)*x(j)+(0.0251e-13)*(y(k))^2+(-54.751e-13)*(x(j))^2+(-14.418e-13)*y(k)*x(j);
%Calcul de LTE à T=800K (pour B)
LTE(j,k)=a0(j,k)+a1(j,k)*800+a2(j,k)*(800^2)+a3(j,k)*(800^3);
%Calcul du coefficient volumétrique d'expansion thermique alpha (pour B)
alpha(j,k)=3*((LTE(j,k))/(800-300));
%Calcul de la densité théorique (pour B)
MmoyMOXkG(j,k)=((1-y(k))/3)*(Mu*10^(-3))+(y(k)/3)*(Mpu*10^(-3))+((2-x(j))/3)*(Mo*10^(-3)); %en kg/mol
MmoyMOXamu(j,k)=((1-y(k))/3)*(Mu)+(y(k)/3)*(Mpu)+((2-x(j))/3)*(Mo); %en g/mol
rhoTh(j,k)=(4*(3-x(j))*(MmoyMOXkG(j,k)))/((a(j,k)^3)*Na);
%Calcul du Bulk Modulus (Pa) (pour B)
BulkModulus(j,k)=((rhoTh(j,k))/3)*(3*(vl(j,k)^2)-4*(vt(j,k))^2);
%Calcul du Cv (pour B)
Cv(j)=3*(3-x(j))*R;
%Calcul du volume molaire (pour B)
Vm(j,k)=((a(j,k)^3)*Na)/(4);
%Calcul du paramètre de Gruneisen (pour B)
Grun(j,k)=(alpha(j,k)*Vm(j,k)*BulkModulus(j,k))/(Cv(j));
%Calcul de la constante B' (pour B)
constB(j,k)=(2.43e-8)/(1-(0.514/Grun(j,k))+(0.228/Grun(j,k)^2));
%Calcul de B
B(j,k)=((1/constB(j,k))*((Grun(j,k)^2)*(4*(3-x(j)))^(2/3))/(MmoyMOXamu(j,k)*(TempDeb(j,k)^3)*VaAngstr(j,k)^(1/3)))*(10^(-2));
%Calcul de la conductivité phonique
CondPhonique(i,j,k)=1/(A(j,k)+B(j,k)*T(i));
%---------------------------------------------------------------------------
%Calcul de la constante d'équlibre de la réaction de disproportiatiation
K(i)=exp((-H/(k1*T(i)))+1);
%Calcul de la fréquence moyenne des phonons (s-1)
v(j,k)=((2*pi*k2)/h)*TempDeb(j,k);
%Calcul des concentrations en trous/electrons
p(i,j,k)=-x(j)+(K(i)-sqrt(4*K(i)*x(j)*(K(i)*(x(j)-2*y(k)+1)+2*(y(k)-x(j))-1)+4*K(i)*y(k)*((y(k)-1)*(K(i)-1))+(K(i)^2)+4*(x(j)^2)))/(2*(K(i)-1));
n(i,j,k)=p(i,j,k)+2*x(j);
% % %Calcul des conductivité électriques dues à Pu et à U (S/m)
condPu(i,j,k)=(((4*(e^2)*v(j,k))/((T(:,i))*a(j,k)*k2))*((n(i,j,k))/(y(k)))*(1-((n(i,j,k))/(y(k)))))*exp(-EhPu/(k1*(T(:,i))));
condU(i,j,k)=(((4*(e^2)*v(j,k))/((T(:,i))*a(j,k)*k2))*((p(i,j,k))/(1-y(k)))*(1-((p(i,j,k))/(1-y(k)))))*exp(-EhU/(k1*(T(:,i))));
% preExpPu(i,j,k)=(((4*(e^2)*v(j,k))/(a(j,k)*k2))*((n(i,j,k))/(y(k)))*(1-((n(i,j,k))/(y(k)))));
% preExpU(i,j,k)=(((4*(e^2)*v(j,k))/(a(j,k)*k2))*((p(i,j,k))/(1-y(k)))*(1-((p(i,j,k))/(1-y(k)))));
% condPu(i,j,k)=((preExpPu(i,j,k))/(T(:,i)))*exp(-EhPu/(k1*(T(:,i))));
% condU(i,j,k)=((preExpU(i,j,k))/(T(:,i)))*exp(-EhU/(k1*(T(:,i))));
% % %Calcul de la conductivité électrique totale (S/m)
condTot(i,j,k)=condPu(i,j,k)+condU(i,j,k);
% condTot(i,j,k)=(((4*(e^2)*v(j,k))/((T(:,i))*a(j,k)*k2))*(n(i,j,k)+p(i,j,k))*(1-n(i,j,k)-p(i,j,k)))*(exp(-0.1/(k1*(T(:,i)))));
% Preexp=((4*(e^2)*v(j,k))/((T(:,i))*a(j,k)*k2))*(n(i,j,k)+p(i,j,k))*(1-n(i,j,k)-p(i,j,k))
%Calcul de la conductivité électronique (W/m/K)
condElectron(i,j,k)=((k2/e)^2)*(T(:,i))*(condPu(i,j)*condU(i,j)/condTot(i,j))*((H/(k1*T(:,i)))^2);
end
end
end

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!