# curve fit a custom polynomial

12 views (last 30 days)
Randy Chen on 26 Oct 2021
Commented: Star Strider on 27 Oct 2021
I have the following 2nd order polynomial in the r-z coordinates: Right now I have four sets of coordinats (r,z), how should I do a curve fit such that I can get an expression for Z in terms of r?
z = [-6.41 -12.4 2.143 102];
r = [13.58 15.7636 12.96 46.6];

Star Strider on 27 Oct 2021
One approach —
syms A B C D r z
sf = A*r^2 + B*z + C*r + D == 0
sf = sfiso = isolate(sf, z)
sfiso = zfcn = matlabFunction(rhs(sfiso), 'Vars',{[A B C D],r})
zfcn = function_handle with value:
@(in1,r)-(in1(:,4)+in1(:,3).*r+in1(:,1).*r.^2)./in1(:,2)
z = [-6.41 -12.4 2.143 102];
r = [13.58 15.7636 12.96 46.6];
B0 = rand(1,4);
nlm = fitnlm(r, z, zfcn, B0)
Warning: The model is overparameterized, and model parameters are not identifiable. You will not be able to compute confidence or prediction intervals, and you should use caution in making predictions.
nlm =
Nonlinear regression model: y ~ F(in1,r) Estimated Coefficients: Estimate SE tStat pValue ________ ________ _______ __________ b1 -0.61998 0.070975 -8.7353 0.072563 b2 2.4979 0.87073 2.8688 0.21353 b3 29.339 1.4571 20.135 0.031591 b4 -275.65 0.16252 -1696.1 0.00037534 Number of observations: 4, Error degrees of freedom: 1 Root Mean Squared Error: 3.89 R-Squared: 0.998, Adjusted R-Squared 0.995 F-statistic vs. constant model: 290, p-value = 0.0415
Bv = nlm.Coefficients.Estimate;
pv = nlm.Coefficients.pValue;
Out = table({'A';'B';'C';'D'},Bv,pv, 'VariableNames',{'Parameter','Value','p-Value'})
Out = 4×3 table
Parameter Value p-Value _________ ________ __________ {'A'} -0.61998 0.072563 {'B'} 2.4979 0.21353 {'C'} 29.339 0.031591 {'D'} -275.65 0.00037534
rv = linspace(min(r), max(r));
zv = predict(nlm, rv(:));
figure
plot(r, z, 'pg')
hold on
plot(rv, zv, '-r')
hold off
grid
xlabel('r')
ylabel('z')
legend('Data','Model Fit', 'Location','best') The Warning was thrown because the number of parameters are not less than the number of data pairs.
Experiment to get different results.
.
##### 2 CommentsShowHide 1 older comment
Star Strider on 27 Oct 2021
As always, my pleasure!
The isolate function was introduced in R2017a.
Use solve instead —
syms A B C D r z
sf = A*r^2 + B*z + C*r + D == 0
sf = sfiso = solve(sf, z)
sfiso = zfcn = matlabFunction(sfiso, 'Vars',{[A B C D],r})
zfcn = function_handle with value:
@(in1,r)-(in1(:,4)+in1(:,3).*r+in1(:,1).*r.^2)./in1(:,2)
z = [-6.41 -12.4 2.143 102];
r = [13.58 15.7636 12.96 46.6];
B0 = rand(1,4);
nlm = fitnlm(r, z, zfcn, B0)
Warning: The model is overparameterized, and model parameters are not identifiable. You will not be able to compute confidence or prediction intervals, and you should use caution in making predictions.
nlm =
Nonlinear regression model: y ~ F(in1,r) Estimated Coefficients: Estimate SE tStat pValue ________ ________ _______ __________ b1 -0.725 0.082997 -8.7353 0.072563 b2 2.9211 1.0182 2.8688 0.21353 b3 34.308 1.7039 20.135 0.031591 b4 -322.35 0.19005 -1696.1 0.00037534 Number of observations: 4, Error degrees of freedom: 1 Root Mean Squared Error: 3.89 R-Squared: 0.998, Adjusted R-Squared 0.995 F-statistic vs. constant model: 290, p-value = 0.0415
Bv = nlm.Coefficients.Estimate;
pv = nlm.Coefficients.pValue;
Out = table({'A';'B';'C';'D'},Bv,pv, 'VariableNames',{'Parameter','Value','p-Value'})
Out = 4×3 table
Parameter Value p-Value _________ _______ __________ {'A'} -0.725 0.072563 {'B'} 2.9211 0.21353 {'C'} 34.308 0.031591 {'D'} -322.35 0.00037534
rv = linspace(min(r), max(r));
zv = predict(nlm, rv(:));
figure
plot(r, z, 'pg')
hold on
plot(rv, zv, '-r')
hold off
grid
xlabel('r')
ylabel('z')
legend('Data','Model Fit', 'Location','best') I like isolate because of the output format, and some of its other characteristics.
.

### More Answers (1)

Rik on 27 Oct 2021
You have two options: rewrite your equation to be a pure quadratic and use polyfit, or use a function like fit or fminsearch on this shape.
If you have trouble implementing either of these two, feel free to comment with what you tried.
##### 2 CommentsShowHide 1 older comment
Rik on 27 Oct 2021
Polyfit will fit a pure polynomial of the form
f(x)=p(1)*x^n +p(2)*x^(n-1) ... +p(n)*x +p(n+1)
That means you can determine the values of -A/B, -C/B, and -D/B with polyfit.
As you may conclude from this: there is no unique solution for your setup, unless you have other restrictions to the values you haven't told yet.
z = [-6.41 -12.4 2.143 102];
r = [13.58 15.7636 12.96 46.6];
p=polyfit(r,z,2)
p = 1×3
0.2482 -11.7452 110.3524