CODEWITHSUNDEEP
matlabcurve-fitting
Ief
Ief

Reputation: 33

Error in lsqcurvefit, cannot calculate correct Jacobian for Arrhenius fit

I am trying to fit an equation of the Arrhenius form to some datapoints using lsqcurvefit. 

    D = D0 * exp( -Ea / ( R * T ));   % Arrhenius equation for curve fitting

D0 and Ea are the values I am looking for. T is temperature and represents X, D is a coefficient and represents Y and R is the Gas constant. Since matlab wouldn't find a solution without supplying a jacobian I calculated the jacobian and included a function, as demonstrated by @m7913d in an earlier post by me ([Matlab curve-fitting won't work for small values (1e-12), what can I do? thanks again!).

When I try to run the code, Matlab returns an Error stating that the supplied Jacobian has the wrong dimensions and that it should have the size of 5-by-2. 

    Error using lsqncommon (line 45)
    Supplied Jacobian is not the correct size:
    the Jacobian matrix should be 5-by-2.

But the Jacobian that corresponds to the main fit-equation is returned by Matlab as a 1-by-4 Matrix. I calculated it in the following way:

    syms D0 Ea R T
    F = D0 * exp(-Ea./(R.* T));
    J = jacobian(F)

    J = [ exp(-Ea/(R*T)), -(D0*exp(-Ea/(R*T)))/(R*T), ... 
    (D0*Ea*exp(-Ea/(R*T)))/(R^2*T), (D0*Ea*exp(-Ea/(R*T)))/(R*T^2)]

But Matlab wouldn't accept this jacobian Matrix to execute the lsqcurvefit operation.

What can I do about it? Did I miss something somewhere? I know that D0 should be something in the order of 1e-5 and Ea is around 170e3.

Any help would be kindly appreciated. Here is a minimum example of the code I use. Please note that this code will result in the Error mentioned above. 

clear all

R1F = [1250 2.5e-11; 1300 2.7e-11; 1350 7.1e-11; 1400 7.2e-11; 1450 1.1e-10];           % test data

 R = 8.3144598;                  % [(kg*m^2) / (s^2 * mol * K)]
 xdata = [R1F(:,1)+272]';
 ydata = R1F(:,2)';

 D0 = 0.1;   % start guess
 Ea = 0.1;   % start guess


 options = optimoptions('lsqcurvefit', 'StepTolerance', 1e-12, ...
'OptimalityTolerance', 1e-12, 'FunctionTolerance', 1e-12, ...
'FiniteDifferenceType', 'central', 'SpecifyObjectiveGradient', true);           
 [X, resnorm, residual, EXITFLAG, OUTPUT] = lsqcurvefit(@(x, xdata) ...
 z(x(1), x(2), xdata, R),[D0 Ea], xdata, ydata, [], [], options);

D0 = X(1);
Ea = X(2);


semilogy(10000./xdata,ydata, '*')
hold on
semilogy(10000./xdata, z(D0, Ea, xdata, R))
hold off


function [F, J] = z(D0, Ea, T, R)
  F = D0 * exp(-Ea./(R.* T));              % function to fit to the datapoints
  J = [ exp(-Ea./(R.*T)), -(D0.*exp(-Ea./(R.*T)))./(R.*T), ...
(D0*Ea*exp(-Ea./(R.*T)))./(R.^2*T), (D0*Ea*exp(-Ea./(R.*T)))./(R.*T.^2)];    % Jacobian of the fit function
end

Upvotes: 0

Views: 490

Answers (2)

Ief
Ief

Reputation: 33

I didn't solve the exact problem I was asking for, but I found a way to work around it. The jacobian I calculated was correct but Matlab just wouldn't find a proper solution for the fit.

So I translated the equation into linear form in log-space. The form of the equation now is:

    log(D) = - (Ea / (ln(10)*R*T) + log(D0)

I was able to fit the parameters using "fit" and calculate values for Ea and D0 from the resulting fit parameters.

So, for me this problem is solved, but if somebody else knows a solution to the initial problem of matlab not finding a solution for the exponential equation, please feel free to share.

Upvotes: 0

m7913d
m7913d

Reputation: 11064

The jacobian should indeed have a size of 5x2:

  • one row for each sample point (length of xdata)
  • one column for each variable you want to fit

Therefore, T should be a column vector and while calculating the jacobian, you should specify to which variables he has to calculate the derivatives. Note that even the order matters!

Upvotes: 1

Related Questions