Solve Non Linear Equations using fsolve

6 views (last 30 days)
Fabian
Fabian on 11 Sep 2012
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.
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
  1 Comment
Star Strider
Star Strider on 11 Sep 2012
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.

Sign in to comment.

Answers (1)

Seth DeLand
Seth DeLand on 11 Sep 2012
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).
  1 Comment
Fabian
Fabian on 11 Sep 2012
Thanks, that actually works.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!