error lsqcurvefit for integrate function

function theta = Test3
% Load data (assuming datapylo.txt exists in your working directory)
data = load('datapylo.txt');
xdata = data(:, 1);
ydata = data(:, 2);

% Initial guess for parameters
k00 = 6^11; % 1/min
E00 = 150000; % J/mol
sig0 = 3000;  % J/mol
theta0 = [k00; E00; sig0];

% Define the model function (vectorized for efficiency)
function y = model(theta, x)
  k0 = theta(1);
  E0 = theta(2);
  sig = theta(3);
  R = 8.314; % J/mol-K  % Constant value, moved outside for clarity
  T = 50;  % Constant temperature value, moved outside for clarity
  pi = 22/7;  % Constant value, consider using symbolic toolbox for pi
  beta = 10;  % Constant value, moved outside for clarity
  
  A = @(E) exp(-E0./(R*T));  % Anonymous function for pre-exponential factor
  f = @(E) 1/(sig*sqrt(2*pi))*exp(-(E-E0).^2. /(2*sig.^2)); % Integrand for distribution
  
  % Perform nested integration numerically (replace with more robust integration if needed)
  C1 =(k0/beta)*(-E0./(R*T) - integral(A,0,T)*k0./beta).*f(x);
  y = integral(C1,0,10);  % Integrate over energy range
end

% Perform curve fitting using lsqcurvefit
[theta_fit, resnorm, residuals] = lsqcurvefit(@model, theta0, xdata, ydata);

% Evaluate the fitted model and plot results
xplot = linspace(min(xdata), max(xdata));
yfit = model(theta0, xplot);

figure;
plot(xdata, ydata, 'o', xplot, yfit);
xlabel('Independent Variable');
ylabel('Dependent Variable');
legend('Data', 'Fitted Curve');

% Display fitting results
disp('Fitted parameters:');
disp(theta_fit);
disp('Residual norm:');
disp(resnorm);
end

this error Failure in initial user-supplied objective function evaluation. LSQCURVEFIT cannot continue.

Thank you for your help

[Torsten]

Which variable in your mathematical description represents xdata, which variable represents ydata ?

[Suravit Naksusuk]

yes

[Torsten]

It was a question ...

Is T the independent variable and alpha(T)/dT the dependent variable (thus T = data(:,1) and dalphadT = data(:,2)) ?

And your mathematical description for dalpha/dT is very sloppy - I cannot figure out for certain where in your formula for dalpha/dT, the variable T is a formal integration variable and where T is a given value.

[Suravit Naksusuk]

Sorry for the not-clear answer.

T and E are dependent variables for d (alpha)/dT.

The xdata is T, and the ydata is d(alpha)/dT.

So I need to integrate dT inside and then find the integration outside for dE.

then compare between d(alpha)/dT,cal and d(alpha)/dT,exp

I attached the file data.

Answer 1

[Torsten]

Replace

% Define the model function (vectorized for efficiency)
function y = model(theta, x)
  k0 = theta(1);
  E0 = theta(2);
  sig = theta(3);
  R = 8.314; % J/mol-K  % Constant value, moved outside for clarity
  T = 50;  % Constant temperature value, moved outside for clarity
  pi = 22/7;  % Constant value, consider using symbolic toolbox for pi
  beta = 10;  % Constant value, moved outside for clarity
  
  A = @(E) exp(-E0./(R*T));  % Anonymous function for pre-exponential factor
  f = @(E) 1/(sig*sqrt(2*pi))*exp(-(E-E0).^2. /(2*sig.^2)); % Integrand for distribution
  
  % Perform nested integration numerically (replace with more robust integration if needed)
  C1 =(k0/beta)*(-E0./(R*T) - integral(A,0,T)*k0./beta).*f(x);
  y = integral(C1,0,10);  % Integrate over energy range
end

by

% Define the model function (vectorized for efficiency)
function y = model(theta, x)
  k0 = theta(1);
  E0 = theta(2);
  sig = theta(3);
  R = 8.314; % J/mol-K  % Constant value, moved outside for clarity
  T = 50;  % Constant temperature value, moved outside for clarity
  pi = 22/7;  % Constant value, consider using symbolic toolbox for pi
  beta = 10;  % Constant value, moved outside for clarity
  for i = 1:numel(x)
    integral_inner = @(E)-k0/beta*integral(@(t)exp(-E./(R*t)),eps,x(i));
    integral_outer = @(E)exp(-E/(R*x(i))).*exp(integral_inner(E)).*exp(-(E-E0).^2/(2*sig^2));
    y(i) = k0/beta/(sqrt(2*pi)*sig)*integral(@(E)integral_outer(E),0,Inf,'ArrayValued',1)
  end
end

[Torsten]

Test it on your computer. The code takes too long to be run here.

I think T must be defined in K instead of degreeC. Further check whether the computation of dalpha/dT is correct. As said, I'm not sure about your mix of formal paramter and value for the variable T and the limit of integration (starting at T=0 with a division by zero in the exp-term).

theta = Test3()

function theta = Test3
% Load data (assuming datapylo.txt exists in your working directory)
data = load('datapylo.txt');
xdata = data(:, 1)+273.15;
ydata = data(:, 2);
% Initial guess for parameters
k00 = 6^11; % 1/min
E00 = 150000; % J/mol
sig0 = 3000;  % J/mol
theta0 = [k00; E00; sig0];
% Perform curve fitting using lsqcurvefit
[theta_fit, resnorm, residuals] = lsqcurvefit(@model, theta0, xdata, ydata,[],[],optimoptions(@lsqcurvefit,'Display','iter','OptimalityTolerance',1e-12))
% Evaluate the fitted model and plot results
xplot = linspace(min(xdata), max(xdata));
yfit = model(theta_fit, xplot);
figure;
plot(xdata, ydata, 'o', xplot, yfit);
xlabel('Independent Variable');
ylabel('Dependent Variable');
legend('Data', 'Fitted Curve');
% Display fitting results
disp('Fitted parameters:');
disp(theta_fit);
disp('Residual norm:');
disp(resnorm);
end


% Define the model function (vectorized for efficiency)
function y = model(theta, x)
  k0 = theta(1);
  E0 = theta(2);
  sig = theta(3);
  R = 8.314; % J/mol-K  % Constant value, moved outside for clarity
  beta = 10;  % Constant value, moved outside for clarity
  y = zeros(size(x));
  for i = 1:numel(x)
    integral_inner = @(E)-k0/beta*integral(@(t)exp(-E./(R*t)),eps,x(i));
    integral_outer = @(E)exp(-E/(R*x(i))).*exp(integral_inner(E)).*exp(-(E-E0).^2/(2*sig^2));
    y(i) = k0/beta/(sqrt(2*pi)*sig)*integral(@(E)integral_outer(E),0,Inf,'ArrayValued',1);
  end
end