curve fit a custom polynomial

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

Answer 1

[Star Strider]

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.

.

[Randy Chen]

Thank you! I was trying out the code myself, but this error occured:

syms A B C D r z ;
sf = A*r^2 + B*z + C*r + D == 0;

sfiso = isolate(sf,z);

zfcn = matlabFunction(rhs(sfiso), 'Vars',{[A B C D],r});

z = [12.96 46.6 13.5 46.5188]
r = [2.143 102 2.41814 101.5];
B0 = rand(1,4);
nlm = fitnlm(r, z, zfcn, B0);

Bv = nlm.Coefficients.Estimate;
pv = nlm.Coefficients.pValue;
Out = table({'A';'B';'C';'D'},Bv,pv, 'VariableNames',{'Parameter','Value','p-Value'});

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


Unrecognized function or variable 'charcmd'.

Error in sym/isolate (line 108)
    throw(CaughtMException);

Error in hw10 (line 4)
sfiso = isolate(sf,z);

[Star Strider]

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.

.

Answer 2

[Rik]

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.

[Randy Chen]

I'm not sure how to rewrite that expression to be pure quadratic. Is it Bz = -Ar^2-Cr-D? But why do I have to rewrite it before using polyfit?

[Rik]

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