fitting implicit non-linear equation

1 Ansicht (letzte 30 Tage)
MATT
MATT am 25 Nov. 2014
Kommentiert: Torsten am 26 Nov. 2014
Dear ALL, I have the data points below, defined by the variables X1,X2,X3,X4:
X1 X2 X3 X4
3.5E-4 0 0.99965 1
4.08935E-4 0.00138 0.99821 0.85588
0.00114 0.03473 0.96413 0.30777
6.43932E-4 0.00989 0.98947 0.54354
0.00145 0.04977 0.94878 0.24193
0.0019 0.08999 0.90811 0.18414
0.00268 0.14986 0.84745 0.13042
0.00268 0.14998 0.84734 0.13056
0.00465 0.2488 0.74655 0.07521
I need to fit these data points using a function of the type
a*b*ln(X4)=W1*X2*(1-X1)+W2*X3*(1-X1)-W3*X2*X3+W4*X2*X3*(1-2*X1)
where a,b are constants (8.31 and 973), and W1,W2,W3,W4 are adjustable parameters. Could you indicate me how to do that? Thank you in advance, Matthieu

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (2)

Torsten
Torsten am 25 Nov. 2014
Use lsqnonlin with f(i) defined as
f(i) = a*b*ln(X4(i))-(W1*X2(i)*(1-X1(i))+W2*X3*(1-X1(i))-W3*X2(i)*X3(i)+W4*X2(i)*X3(i)*(1-2*X1(i)))
Best wishes
Torsten.
  2 Kommentare
MATT
MATT am 25 Nov. 2014
Bearbeitet: MATT am 25 Nov. 2014
Hello, and thank. I am not sure I understand how to input the indices:
I define my function as:
function F = FIT(W, X1, X2, X3, X4)
i=1:9;
F = 8.31*973*ln(X4(i))-(W(1)*X2(i)*(1-X1(i))+W(2)*X3(i)*(1-X1(i))-W(3)*X2(i)*X3(i)+W(4)*X2(i)*X3(i)*(1-2*X1(i)));
end
Then I launch the following command:
W0=[0.10,1,40,1];
[W] = lsqnonlin(@FIT,W0, X1, X2, X3, X4);
This does not seem to work, I guess it inputs all of the X(i) rather than fitting the data. Any advice on how to debbug it? Best, Matthieu
Torsten
Torsten am 26 Nov. 2014
As I see now, your function is linear in the unknown Parameters.
So you can determine W simply by
A=zeros(9,4);
B=zeros(9);
for ii=1:9
A(ii,1)=X2(ii)*(1-X1(ii));
A(ii,2)=X3(ii)*(1-X1(ii));
A(ii,3)=-X2(ii)*X3(ii);
A(ii,4)=X2(ii)*X3(ii)*(1-2*X1(ii));
B(ii)=a*b*log(X4(ii));
end
W=A\B;
Best wishes
Torsten.

Melden Sie sich an, um zu kommentieren.

arich82
arich82 am 26 Nov. 2014
I question the physicality of your fitting equation:
Since you have a log on the left-hand side, and your data includes X4(1) == 1 and X2(1) == 0 exactly, it seems like W2 should necessarily be 0 for the equation to remain valid...
In your data, it seems like X3 is very nearly (1 - X2). If you use this substitution, you can plot your data (and your fit) with plot3.
Feel free to play with the code below. I used the backslash operator for fitting, which I believe results in a least-squares fit via QR (but don't hold me to that):
First, some preliminary definitions.
%%%%
%%data
%%%%
clear; close all; clc;
data = [...
0.00035, 0, 0.99965, 1; ...
0.000408935, 0.00138, 0.99821, 0.85588; ... % swapped rows
0.000643932, 0.00989, 0.98947, 0.54354; ... % swapped rows
0.00114, 0.03473, 0.96413, 0.30777; ...
0.00145, 0.04977, 0.94878, 0.24193; ...
0.0019, 0.08999, 0.90811, 0.18414; ...
0.00268, 0.14986, 0.84745, 0.13042; ...
0.00268, 0.14998, 0.84734, 0.13056; ...
0.00465, 0.2488, 0.74655, 0.07521; ...
];
% Note: swapping rows 3 & 4 dowsn't affect the fitting, but makes
% ploting lines cleaner
X1 = data(:, 1);
X2 = data(:, 2);
X3 = data(:, 3);
X4 = data(:, 4);
a = 8.31;
b = 973;
% re-parameterize (unnecessary, but might simplify some later algebra)
T = a*b*log(X4);
X = X2;
Y = 1 - X1;
Z = X3;
% Note: Z = X3 ~~ (1 - X2) == (1 - X)
For the first fit, we'll assume X3 = 1 - X1
%%%%
%%w: fit assuming X3 = (1 - X2), using all four w's
%%%%
% T = w2*Y + (w1 - w2 + 2*w4)*X.*Y + (w3 + w4)*(X.^2 - X) - 2*w4*X.^2.*Y;
M = [Y, X.*Y, X.^2 - X, X.^2.*Y];
m = M\T;
Tm = M*m;
X4m = exp(Tm/a/b);
w(4) = m(4)/2;
w(3) = m(3) - w(4);
w(2) = m(1);
w(1) = m(2) + w(2) - 2*w(4);
w = w(:);
err_m = (X4 - X4m)./X4m;
w
err_m
Now let's see what happens if we enforce the (perhaps more physical) constraint W2 = 0
%%%%
%%w2: fit assuming X3 = (1 - X2), using w2 = 0
%%%%
% T = (w1 + 2*w4)*X.*Y + (w3 + w4)*(X.^2 - X) - 2*w4*X.^2.*Y;
M = [X.*Y, X.^2 - X, X.^2.*Y];
m2 = M\T;
Tm2 = M*m2;
X4m2 = exp(Tm2/a/b);
w2(4) = m(3)/2;
w2(3) = m(2) - w2(4);
w2(2) = 0;
w2(1) = m(2) - 2*w2(4);
w2 = w2(:);
err_m2 = (X4 - X4m2)./X4m2;
w2
err_m2
Now lets use the full equation (no assumption on X3)
%%%%
%%W: fit using all four w's (no assumptions)
%%%%
%
% TW = ...
% W(1)*X2.*(1 - X1) ...
% + W(2)*X3.*(1 - X1) ...
% - W(3)*X2.*X3 ...
% + W(4)*X2.*X3.*(1 - 2*X1);
M = [X2.*(1 - X1), X3.*(1 - X1), -X2.*X3, X2.*X3.*(1 - 2*X1)];
W = M\(a*b*log(X4));
TW = M*W;
X4W = exp(TW/a/b);
err_W = (X4 - X4W)./X4W;
W
err_W
And finally, using W2 == 0 (no assumption on X3)
%%%%
%%W: fit using w2 = 0 (no additional assumptions)
%%%%
M = [X2.*(1 - X1), -X2.*X3, X2.*X3.*(1 - 2*X1)];
W2 = M\(a*b*log(X4));
TW2 = M*W2;
X4W2 = exp(TW2/a/b);
err_W2 = (X4 - X4W2)./X4W2;
W2 = [W2(1), 0, W2(2), W2(3)].';
W2
err_W2
results
%%%%
%%compare results
%%%%
compare_W = [w, w2, W, W2];
compare_W
compare_err = 100*[err_m, err_m2, err_W, err_W2];
compare_err
hf = figure('WindowStyle', 'docked');
ha = axes;
hp = plot3(X1, X2, X4, X1, X2, X4m, X1, X2, X4m2, X1, X2, X4W, X1, X2, X4W2)
xlabel('X1')
ylabel('X2')
zlabel('X4')
legend('X4', 'X4m', 'X4m2', 'X4W', 'X4W2')
grid on;
title('compare fits')
You can use your judgment to see if this %error is reason able over the region of interest.
>> compare_W =
1.0e+08 *
3.813510730697298 -0.000491562357117 0.025845576531701 0.032553401777266
-0.000018277710575 0 -0.000018555355233 0
2.867011340403814 0.959509492664202 1.017723304845835 1.273754444212918
-0.947009230361176 0.960001055021319 0.991441774206119 1.240781494391939
>> compare_err =
25.3541 0 25.7751 0
9.5807 -12.2780 9.8547 -12.3620
-18.5631 -32.9458 -18.6725 -33.2918
-19.8587 -26.6898 -20.2522 -27.2932
-5.6106 -6.3826 -5.7995 -6.7560
19.3484 26.9156 19.2132 26.7386
0.7148 2.1690 0.9361 2.5826
0.2742 1.5628 0.5416 2.0342
-1.9625 -3.4741 -2.0987 -3.7206

Kategorien

Mehr zu Interpolation finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by