Curve fitting with a constrained y value to Zero
8 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Khaled Bahnasy
le 13 Août 2020
Modifié(e) : Matt J
le 13 Août 2020
I need to fit a curve ( Second-order polynomial ) with a constraint that a specific y-value equal to Zero
4.4 2.367224698
21.1 0
37.8 -1.857318083
54.4 -3.276015126
X & Y values as an example attached X = [ 4.4 21.1 37.8 54.4 ]
I want to fit the cuve where the y-value at x= 21.1 equal to zero
I am new to matlab and i have tried Curve fitting toolbox, I think it is not provided as a constraint in the toolbox
0 commentaires
Réponse acceptée
John D'Errico
le 13 Août 2020
Modifié(e) : John D'Errico
le 13 Août 2020
This is quite easy, actually.
xy = [4.4 2.367224698
21.1 0
37.8 -1.857318083
54.4 -3.276015126];
x = xy(:,1);
y = xy(:,2);
Now, you want to force a quadratic polynomial to go through y==0, at x == 21.1.
The simple solution is to fir a quadratic model of the form
y = a1*(x - 21.1) + a2*(x - 21.1)^2
So with effectively no constant term. What happens when x == 21.1? WE GET ZERO! See how nicely it works? Yes, it is true, that IF you expand it out there would be a constant term. And that is why you use this form for the model.
That quadratic model is not difficult to fit, either.
a12 = [x - 21.1,(x-21.1).^2]\y;
a1 = a12(1)
a1 =
-0.126909925659631
a2 = a12(2)
a2 =
0.000863233648946335
The model is simple to evaluate.
ypred = a1*(x - 21.1) + a2*(x-21.1).^2
ypred =
2.36014299087048
0
-1.8786485261612
-3.26886936348561
As you can see, it passes EXACTLY through the point (21.1,0).
Could we have gotten this in the form of a 2 parameter model, so with 3 coefficients? Well yes. That would take slightly more effort, but still quit doable.
Perhaps simplest, if we wanted to find the 3 parameter quadratic that has the indicated behavior from a1 and a2, we could just expand the polynomial. Pencil and paper will suffice. Thus is we have this function:
a1*(x - 21.1) + a2*(x-21.1).^2
then in the form
b0 + b1*x + b2*x.^2
We would have
b2 = a2
b1 = a1 - 42.2*a2
b0 = 445.21*a2 - 21.1*a1
Again, if we use this to predict the value for the vector x, we see:
b0 + b1*x + b2*x.^2
ans =
2.36014299087048
-6.10622663543836e-16
-1.8786485261612
-3.26886936348562
The second element of that vector is non-zero, but only by an amount that corresponds to floating point crap in the last significant bits of the number. 6e-16 is effectively zero here.
3 commentaires
John D'Errico
le 13 Août 2020
Funny. As you were asking me how to write that in terms of 3 coefficients, I was explaining how to do that by amending my answer.
John D'Errico
le 13 Août 2020
Note that if I am feeling too lazy to do the pencil and paper (not uncommon for me) then I might have let MATLAB do the work.
syms a1 a2 X
vpa(expand(a1*(X - 21.1) + a2*(X-21.1).^2))
ans =
445.21*a2 - 21.1*a1 + X*a1 - 42.2*X*a2 + X^2*a2
Now by collecting the terms, we should get the same conversion I wrote.
Plus de réponses (2)
Serhii Tetora
le 13 Août 2020
Modifié(e) : Serhii Tetora
le 13 Août 2020
clear;close;clc
x = [4.4 21.1 37.8 54.4 ];
y = [2.367224698 0 -1.857318083 -3.276015126];
w = [1 1000 1 1];
[xData, yData, weights] = prepareCurveData( x, y, w );
% Set up fittype and options.
ft = fittype( 'poly2' );
opts = fitoptions( 'Method', 'LinearLeastSquares' );
opts.Weights = weights;
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
% Plot fit with data.
h = plot( fitresult, xData, yData, 'o' );
legend( h, 'y vs. x', 'fit', 'Location', 'NorthEast');
% Label axes
xlabel('x');
ylabel('y');
grid on
Matt J
le 13 Août 2020
Modifié(e) : Matt J
le 13 Août 2020
Using lsqlin,
x0=21.1;
x = [4.4 37.8 54.4 ].';
y = [2.367224698 -1.857318083 -3.276015126].';
p=lsqlin(x.^[2,1,0],y,[],[],x0.^[2,1,0],0)
Check:
>> polyval(p,21.1)
ans =
-4.4409e-16
2 commentaires
Matt J
le 13 Août 2020
Modifié(e) : Matt J
le 13 Août 2020
The approximation error given by my approach is at the limit of floating point precision. It's meaningless to aspire beyond that. The value 21.1 doesn't even have an exact representation in a binary floating point computer. To 47 decimal places, the number that the computer is really holding when you enter x0=21.1 is,
'21.10000000000000142108547152020037174224853515625'
Voir également
Catégories
En savoir plus sur Get Started with Curve Fitting Toolbox dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!