MATLAB nlinfit with multiple unknowns

Question

I am trying to fit a curve using nlinfit, however, I have 3 unknowns.

The equation that I am trying to fit is the laminar velocity profile in a tube.

u = umax( 1 - ( r / R ) ^ 2 ), where u is the velocity, umax is the centerline velocity, r is the distance from the centerline and R is the radius of the tube.

The problem is that umax, r and R are unknowns.

http://www.mne.psu.edu/cimbala/learning/fluid/CV_Momentum/pipe_eg.gif

From the image above, the centerline is the axis. I have included code for 1 set of data:

clc
clear all

XMean = [0.13518 
0.1599
0.17035
0.18053
0.18849
0.19577
0.19373
0.18781
0.17245
0.15728
0.13404
0.10981];

r = [-5.5000
   -4.5000
   -3.5000
   -2.5000
   -1.5000
   -0.5000
    0.5000
    1.5000
    2.5000
    3.5000
    4.5000
    5.5000];

plot( XMean, r )

I am not sure how to format the equation handle:

Eqn = @(u,y) u(1).*( 1 - ( ( u(2) - y ) / u(3) ).^2 );

[ beta, R, J, CovB, MSE ] = nlinfit( YMean, r, Eqn, Alpha );

I would appreciate any help trying to fit this data. Thanks for your time.

Answer 1

In order to compute the rmse it is neccessary to have a maximum of 1 y value per x value, mathematically speaking a "function". You can plot later again the way your sketch shows it. But for the computation: rotate by 90 degree.

For the optimization, you need an initiall guess (your Alpha), which is not defined. However, with some minor changes in the code below, one can get the optimal parameters. Also it might be better readable. The plot is shown below.

plot(r, XMean ,'ko', 'markerfacecolor', 'black') % // need a "function" (max 1y value for 1x-value)
Eqn = @(par, x) par(1).*( 1 - ( (par(2) - x) / par(3) ).^2 );
par0 = [1, 1, 1]; % // an initial guess as start for optimization
par_fit = nlinfit(r, XMean, Eqn, par0 ); % // optimize
hold on
r_fit = linspace(min(r), max(r)); % // more points for smoother curve
plot(r_fit, Eqn(par_fit, r_fit), 'r')

The optimal parameters i get returned are

par_fit =

    0.1940   -0.4826    9.0916

You might want to rotate the plot again to get it into the right orientation, as in your sketch.

Additional Info: polyfit() instead of nlinfit

Since one can compute this flow profile analytically, one knows in advance that is is a parabola. You can also use polyfit to fit a parabola.

par_poly_fit = polyfit(r, XMean, 2)

% // it yields the values:
par_poly_fit =

-0.0023   -0.0023    0.1934

Plotting that paraobla gives a line pretty much identical to the line, which we got with the optimizer, but it is just more stable since it does not depend on an initial value.

 plot(r_fit, polyval(par_poly_fit, r_fit),'g--'