首页 文章

lsqcurvefit中的错误，无法为Arrhenius拟合计算正确的雅可比行列式

提问于
浏览 次
0

我试图使用lsqcurvefit将Arrhenius形式的方程拟合到某些数据点 .

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

D0和Ea是我要找的值 . T是温度，表示X，D是系数，表示Y，R是气体常数 . 由于matlab在没有提供jacobian的情况下找不到解，我计算了jacobian并包含了一个函数，正如@ m7913d在我之前的帖子中所证明的那样（[Matlab curve-fitting won't work for small values (1e-12), what can I do?再次感谢！） .

当我尝试运行代码时，Matlab会返回一个错误，指出所提供的雅可比行星具有错误的尺寸，并且其大小应为5乘2 .

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

但是对应于主拟合方程的雅可比行列式由Matlab作为1乘4矩阵返回 . 我用以下方式计算它：

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

但是Matlab不会接受这个雅可比矩阵来执行lsqcurvefit操作 .

我能做些什么呢？我在某处错过了什么吗？我知道D0应该是1e-5的数量级，Ea是170e3左右 .

任何帮助将不胜感激 . 这是我使用的代码的最小示例 . 请注意，此代码将导致上述错误 .

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
matlab curve-fitting

2 回答

  • 1

    jacobian的大小确实应该是 5x2 ：

    • 每个样本点一行（长度 xdata ）

    • 您希望适合的每个变量都有一列

    因此， T 应该是列向量，在计算雅可比时，您应该指定他有哪些变量来计算导数 . 请注意，即使订单也很重要！

    回复于
  • 0

    我没有解决我要求的确切问题，但我找到了解决问题的方法 . 我计算的jacobian是正确的，但Matlab只是找不到合适的解决方案 .

    所以我将方程式转换为对数空间中的线性形式 . 等式的形式现在是：

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

    我能够使用“拟合”拟合参数，并从得到的拟合参数计算Ea和D0的值 .

    所以，对我来说这个问题已经解决了，但是如果其他人知道matlab最初问题的解决方案没有找到指数方程的解决方案，请随意分享 .

    回复于

相关问题