2017-09-18 4 views
0

J'essaie d'ajuster une équation de la forme Arrhenius à certains points de données en utilisant lsqcurvefit.Erreur dans lsqcurvefit, ne peut pas calculer jacobien correct pour Arrhenius fit

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

D0 et Ea sont les valeurs que je recherche. T est la température et représente X, D est un coefficient et représente Y et R est la constante de gaz. Puisque matlab ne trouverait pas de solution sans fournir un jacobien, j'ai calculé le jacobien et inclus une fonction, comme démontré par @ m7913d dans un post précédent par moi ([Matlab curve-fitting won't work for small values (1e-12), what can I do? merci encore!).

Lorsque j'essaie d'exécuter le code, Matlab renvoie une erreur indiquant que le jacobien fourni a les mauvaises dimensions et qu'il doit avoir la taille de 5 par 2.

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

Mais le jacobien qui correspond à la forme principale équation est renvoyée par Matlab en tant que 1-en-4 Matrix. Je calculé de la manière suivante:

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

Mais Matlab n'accepter cette matrice jacobienne pour exécuter l'opération de lsqcurvefit.

Que puis-je faire à ce sujet? Ai-je manqué quelque chose quelque part? Je sais que D0 devrait être quelque chose dans l'ordre de 1e-5 et Ea est autour de 170e3.

Toute aide serait appréciée. Voici un exemple minimal du code que j'utilise. Veuillez noter que ce code entraînera l'erreur mentionnée ci-dessus.

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 

Répondre

1

Le jacobian devrait en effet avoir une taille de 5x2:

  • une ligne pour chaque point d'échantillon (longueur de xdata)
  • une colonne pour chaque variable que vous souhaitez adapter

Par conséquent, T doit être un vecteur de colonne et lors du calcul de la jacobian , vous devez spécifier à quelles variables il doit calculer les dérivés . Notez que même l'ordre compte!

+0

Merci @ m7913d pour votre aide. Je n'étais toujours pas capable de le comprendre même si j'ai réussi à obtenir le jacobien correct (j'espère). Matlab résout pour deux variables mais les valeurs ne sont pas correctes et l'intrigue vérifie que c'est faux.Je ne sais toujours pas si j'ai calculé le jacobien correctement. Cela vous dérangerait-il de le calculer et de le partager ici? Je dois résoudre pour 'Ea' et 'D0'. – Ief

0

Je n'ai pas résolu le problème exact que je demandais, mais j'ai trouvé un moyen de contourner le problème. Le jacobien que j'ai calculé était correct mais Matlab ne trouverait pas de solution appropriée pour l'ajustement. J'ai donc traduit l'équation en forme linéaire dans l'espace logarithmique. La forme de l'équation est maintenant:

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

j'ai pu adapter les paramètres à l'aide « adapter » et de calculer les valeurs de Ea et D0 des paramètres d'ajustement résultant. Donc, pour moi, ce problème est résolu, mais si quelqu'un d'autre connaît une solution au problème initial de matlab ne trouvant pas de solution pour l'équation exponentielle, n'hésitez pas à partager.