Fit with an integral involved

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.

[Steven Lord]

Can you show us the function you wrote and how you called lsqcurvefit? It may be that your code only need small modifications to work, and that would be easier than starting from scratch.

In addition, can you clarify what "does not work" mean in this context?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support using the telephone icon in the upper-right corner of this page so we can investigate.

[Daniele Sonaglioni]

Sorry but i have answered to my own question so you can see my answer below.

Sorry, my mistake!

[Daniele Sonaglioni]

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);

Answer 1

[Star Strider]

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.

[Daniele Sonaglioni]

I'm sorry because i have tried but it does not work with my data.

Below is the error:

Error using lsqcurvefit (line 263)

Function value and YDATA sizes are not equal.

Error in trial (line 8) %i have created a function_m_file for the subroutine you have written for defining the fitting function

B = lsqcurvefit(@(b,T) integralfit(b,T,R,nu,T0), rand(4,1), T, y)

I have simply copied and pasted your code and used with my data.

Before i have posted an answer to my question with an extract of my data and my previous code.

[Star Strider]

This runs without error with the function I posted:

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];

R = 8.314;                                                                  % Use Actual Value
nu = 100;                                                                 % Use Actual Value
T0 = T_tot(1);                                                              % Use Actual Value

B = lsqcurvefit(@(b,T)integralfit(b,T,R,nu,T0), rand(4,1)*1E+3, T_tot, Cp_tot)

figure
plot(T_tot, Cp_tot, 'p')
hold on
plot(T_tot, integralfit(B,T_tot,R,nu,T0), '-r')
hold off
grid

The problem is that it does not provide a reasonable estimate of the parameters, and that is likely due to the initial parameter estimates.

I have been running it using ga to estimate the parameters, however I still cannot get a decent fit. If you have any information on at least what the magnitudes of the relevant parameters should be, please share it so that I can scale the ga population matrix appropriately.

[Daniele Sonaglioni]

I am finding the same problem.

Reffering to your mapping of the parameter, i can give the magnitude of the parameters:

b(1) = E1 = 1e5

b(2) = E2 = 1e5

b(3) = Tstar = 10/100

b(4) = Tdot = 10/100

This is only a rough estimation but they are the values i expect to find from this fitting

[Daniele Sonaglioni]

With "10/100" i mean a value between 10 and 100

[Star Strider]

I still cannot get a decent fit.

I checked my code and it appears to be essentially what you posted, so check it to be certain it does what you want it to and that it reproduces what you posted. The ‘k3’ and ‘K’ functions should have access to the parameters (they woiuld throw errors if they did not), so that is not the problem. I also do not see any errors in the ‘y’ assignment.

Unless I am not understanding something about the equations you posted, my code should produce appropriate parameter estimates with ga, since it searches the entire parameter space for the best parameter set.

If you find any errors in my code, please post back here with appropriate corrections so that I can continue to see if a good set of parameters exists for it.

I am stopping here, since this does not currently appear to be a solvable problem.

[Daniele Sonaglioni]

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]

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.