Solve Non Linear Equations using fsolve

Hi I´m trying to estimate the center and radius of a sphere. I know the location of certain points f the surface of the sphere, so using this information and the general equation of a sphere I can generate a non-linear equation system and solve it to find the center and radius.

Here is the function I´m using

function F = centerEstimator(sphere)

%%Sphere general equation (x-h)^2 + (y-j)^2 + (z-k)^2 - r^2 = 0
%%where h, j and k are the center coordintes and r is the radius of the
%%sphere

cx = sphere(1);
cy = sphere(2);
cz = sphere(3);
r = sphere(4);

global n;
global X;

for i=1:n
    F(i) = sqrt((X(1,i)-cx)^2 + (X(2,i)-cy)^2 + (X(3,i)-cz)^2) - r;
end

Where X are the coordinates of some points on the surface of the sphere and n is the number of points (1600).

Then I try to use fsolve to estimate cx, cy, cz and r. Using the mean values of the positions of the know points as initial values for fsolve.

options=optimset('TolFun' , 1e-6,'Display','iter', 'TolX' , 1e-8, 'MaxFunEvals' , 5000, 'MaxIter' , 5000);
[estimatedSphere r] = fsolve(@centerEstimator,sphere,options);

As an output I get that fsolve does not found a solution:

Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead. > In fsolve at 314 In CalibrationValidate_v3 at 117

                                        First-Order                    Norm of 
 Iteration  Func-count    Residual       optimality      Lambda           step
     0           5          814796        3.93e+05         0.01
     1          10          230434        1.69e+05        0.001        4.05635
     2          15         10156.8        3.54e+04       0.0001        1.33507
     3          20         4464.92             499        1e-05       0.160848
     4          25         4463.72           0.111        1e-06     0.00240457
     5          30         4463.72        3.48e-07        1e-07    5.37669e-07

No solution found.

fsolve stopped because the problem appears regular as measured by the gradient, but the vector of function values is not near zero as measured by the selected value of the function tolerance.

criteria details

I really don´t know what does this means. I´m doing something wrong? Fsolve is not capable of solving this problem? Does the values I get after this output is the best solution for this problem?

Thanks

[Star Strider]

I recently experimented with finding the center of a circle (for fun, actually), and developed this:

    % Generate data: 
    Cc = [2  -3];
    Cr = 5;

    a = [0:0.3:(pi/2)]';
    for k1 = 1:length(a)
        Ccal(k1,:) = Cc + Cr*[cos(a(k1))  sin(a(k1))];
    end

    % Solve circle equation ‘(x - B(1))^2 + (y - B(2))^2 = B(3)^2’ for y: 
    %   B(3)^2 - (x - B(1))^2 = (y - B(2))^2
    %   sqrt( B(3)^2 - (x - B(1))^2 ) = (y - B(2))
    %   B(2) + sqrt( B(3)^2 - (x - B(1))^2 ) = y

    RSc = @(B,x) B(2) + sqrt( (B(3).^2) - (x - B(1)).^2);
    B0 = rand(3,1)*1;
    [B_est] = lsqcurvefit(RSc, B0, Ccal(:,1), Ccal(:,2) )

    B_est0 = sign(real(B_est)).*abs(B_est)
    Y_fit = RSc(B_est, Ccal(:,1));
    Y_fit = sign(real(Y_fit)).*abs(Y_fit);

    figure(1)
    plot(Ccal(:,1), Ccal(:,2), '-b')
    hold on
    plot(Ccal(:,1), Y_fit, '+g')
    plot(B_est0(1), B_est0(2), 'pr')
    hold off
    grid
    axis equal

It works for circles (although like all nonlinear regressions it occasionally gets lost in the wilderness), so it usually requires a few restarts before it converges on appropriate parameter estimates. The sqrt function generates complex numbers but with imaginary parts on the order of 1E-6 times the real parts, so I put in a couple lines of code to correct for that. Since the nonlinear regression functions nlinfit and lsqcurvefit can take matrix independent variables and can fit matrix dependent variables, several variations for fitting sphere data are possible. The only requirement is that the data your objective function returns matches the dimensions of the dependent data you ask the nonlinear regression function to fit.

An analogous solution to fit a sphere would likely work. This is not really an ‘Answer’ so I'm posting it as a comment instead. There are likely many other ways to solve your problem.

Answer 1

[Seth DeLand]

From what I can tell this problem has 1600 equations and only 4 unknowns, which would mean it is overdetermined and likely doesn't have a solution. Perhaps instead of trying to solve so that the sphere passes exactly through all 1600 points (which probably isn't possible), could you try to minimize the error between the sphere and the points?

I think it would be pretty easy in this case to swap out fsolve for something like lsqnonlin which would minimize the sum of the squared error for your problem.

fsolve tries to find the values that make each element of F in your problem 0. lsqnonlin would try to minimize sum(F.^2).

[Fabian]

Thanks, that actually works.