Fit with an integral involved

3 次查看（过去 30 天）
Daniele Sonaglioni
评论： Star Strider 2021-3-3
Hi everybody,
i have some problem in fitting my experimental data with a function involving an integral in it.
My fitting function is the subsequent:
where
and
where E1,E2,T^* and T^' are my fitting parameters, while T is my abscissa.
I have tried defining a unique function in lsqcurvefit but it does not work. I have some difficulties in defining the integral in the fitting function.
Can someone help me with the code?
Thaks a lot.
  3 个评论
显示 1更早的评论
Daniele Sonaglioni
Daniele Sonaglioni 2021-3-3
Sorry but i have answered to my own question so you can see my answer below.
Sorry, my mistake!
Daniele Sonaglioni
Daniele Sonaglioni 2021-3-3
Hello,
the code i have written works, in the sense that Matlab do not sends any error message but the result of the fitting procedure is not what i want.
Here is an extract of my x and y data:
T_tot=[98.2810 99.2830 100.2800 101.2800 102.2800 103.2800 104.2800 105.2800 106.2800 107.2800 108.2700 109.2700 110.2700 111.2700 112.2700 113.2700 114.2700 115.2700 116.2700 117.2700 118.2700 119.2600 120.2600 121.2600 122.2600]
Cp_tot=[0 0.0003 0.0013 0.0017 0.0027 0.0048 0.0064 0.0087 0.0099 0.0138 0.0159 0.0178 0.0196 0.0213 0.0211 0.0205 0.0187 0.0152 0.0117 0.0083 0.0056 0.0024 0.0013 0.0001 -0.0008]
Below i report the code i have written in which i want to fit what i call T_tot and Cp_tot:
clear
load 'prova_100Ks.txt'
prova_100Ks;
T1=prova_100Ks(:,1);
Cp1=prova_100Ks(:,2);
T_before=[50:.2:min(T1)];
le=length(T_before);
Cp_before=zeros(1,le);
T_after=[max(T1):0.2:150];
l=length(T_after);
Cp_after=zeros(1,l);
T_tot=[T_before'; T1; T_after'];
Cp_tot=[Cp_before';Cp1;Cp_after']; %heat flow in mW
Cp_tot=(Cp_tot.*1e-3)./(100*1e-6);%Cp in J/g*k
Cp_tot=(Cp_tot.*14300)./1000;%in J/K*mol
R=8.314;
rate=100;
integral_term = @(a,b,c,d,x) integral(@(x) exp(-a*(1./x-1/b)).*exp(-c*(1./x-1/d)) , min(x), max(x));
fun = @(x,T_tot) (x(3)/rate).* exp(-x(1)*(1./T_tot-1/x(2))).*exp(-x(3)*(1./T_tot-1/x(4))).* exp(-(1/rate)* integral_term(x(1),x(2),x(3),x(4), T_tot));
fun2 = @(x,Xdata) arrayfun(@(T_tot) fun(x,T_tot), Xdata);
x0 = [1,1,1,1];
x = lsqcurvefit(fun2, x0, T_tot, Cp_tot);

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2021-3-3
Try this:
function y = integralfit(b,Tv,R,nu,T0)
% % % MAPPING: b(1) = E1, b(2) = E2, b(3) = Tstar, b(4) = Tdot
for k = 1:numel(Tv)
k3 = @(T) exp(-b(2).*(1./T - 1./b(3))./R);
K = @(T) exp(-b(1).*(1./T - 1./b(4))./R);
y(k) = (b(1)/nu) .* k3(Tv(k)).* K(Tv(k)) .* exp(-(1./nu).*integral(@(T) k3(T).*K(T), T0, Tv(k), 'ArrayValued',1));
end
end
T = 1:10; % Use Actual Vector
y = rand(size(T)); % Use Actual Vector
R = 9; % Use Actual Value
nu = 7; % Use Actual Value
T0 = T(1); % Use Actual Value
B = lsqcurvefit(@(b,T)integralfit(b,T,R,nu,T0), rand(4,1), T, y)
figure
plot(T, y, 'p')
hold on
plot(T, integralfit(B,T,R,nu,T0), '-r')
hold off
grid
It runs without error and appears to produce appropriate results on a random vector and random extra parameters. I obviously cannot test it with your data to see if it is producing appropriate parameter estimates.
It is somewhat slow because of the loop, however there is no way to correct for that since ‘T’ is the upper limit of integration, and that must be a scalar for integral to work correctly.
Remember to save the ‘integralfit’ function as integralfit.m on your MATLAB search path before you use it.
  7 个评论
显示 5更早的评论
Daniele Sonaglioni
Daniele Sonaglioni 2021-3-3
I have improved my fit by the use of fminsearch and now i have obtained something good but the quality of the fitting parameters are not the best at all.
Anyway, i want to really thank you for your help!
Below i show the code i have added, just for completness:
b=[1e5,1e5,100,100];
f=@(b) integralfit(b,T,R,nu,T0);
z=@(b) sum((abs(y-f(b))).^2);
B=fminsearch(z,b);
Star Strider
Star Strider 2021-3-3
As always, my pleasure!
I am glad that you were able to estimate the parameters correctly, since they eluded me, even using the genetic algorithm.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 t Location-Scale Distribution 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by