Fit with an integral involved

6 views (last 30 days)
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 Comments
Daniele Sonaglioni
Daniele Sonaglioni on 3 Mar 2021
Sorry but i have answered to my own question so you can see my answer below.
Sorry, my mistake!
Daniele Sonaglioni
Daniele Sonaglioni on 3 Mar 2021
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);

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 3 Mar 2021
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 Comments
Daniele Sonaglioni
Daniele Sonaglioni on 3 Mar 2021
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 on 3 Mar 2021
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!