fsolve vs. lsqnonlin: Three intersecting spheres

13 views (last 30 days)
I'm trying to determine the coordinates of the intersection point of three sphere surfaces utilizing either of fsolve and lsqnonlin. I chose these two because others have posted their code for comparable problems. The code for both seems to work. However, judging by my simple distance checks between the sphere centers and the solutions, lsqnonlin seems to result in a rather big variance (e.g. distances 11.7,10.0,12.97 for radii of 10) whereas fsolve seems to be spot-on (distances == 10). I'd therefore use fsolve although slight inaccuracies of ~+5% would be ok(-5% wouldn't be ok). On the other hand, fsolve doesn't support the bounds I need (e.g. x(1) > 0; [0;-Inf;-Inf] for lsqnonlin).
Q: Can I make fsolve return both intersection points (if both exist) so I can select the correct one manually? Or can I make lsqnonlin more accurate?
Code I use to call fsolve:
coord = [x1,y1,z1;x2,y2,z2;x3,y3,z3];
r = [r1;r2;r3]
param = %irrelevant for this;
start=[0;0;0];
[x,fval]= fsolve(@(unknown)surfaceFunctions2(unknown, coord, r, param), start);
%Distance check:
dist1=sqrt((x(1)-y1)^2+(x(2)-x1)^2+(x(3)-z1)^2);
dist2=sqrt((x(1)-y2)^2+(x(2)-x2)^2+(x(3)-z2)^2);
dist3=sqrt((x(1)-y3)^2+(x(2)-x3)^2+(x(3)-z3)^2);
surfaceFunctions2 contains:
function F = surfaceFunctions2(unknown, coord, r, opt)
y = unknown(1);
x = unknown(2);
z = unknown(3);
%Surface functions for the spheres
f1 = (x - coord(1,1))^2 + (y - coord(1,2))^2 + (z - coord(1,3))^2 - r(1)^2;
f2 = (x - coord(2,1))^2 + (y - coord(2,2))^2 + (z - coord(2,3))^2 - r(2)^2;
f3 = (x - coord(3,1))^2 + (y - coord(3,2))^2 + (z - coord(3,3))^2 - r(3)^2;
F = [ f1 ; f2 ; f3];
end;
I use this to call lsqnonlin:
coord = [x1,y1,z1;x2,y2,z2;x3,y3,z3];
r = [r1;r2;r3]
param = %irrelevant for this;
start=[0;0;0];
lb = [0; -Inf;-Inf];
ub = [Inf;Inf;Inf];
[x, fval] = lsqnonlin(@(x)sphereSurfaces(x, coord,r,param), start, lb, ub, opt);
%Distance check:
dist1=sqrt((x(1)-y1)^2+(x(2)-x1)^2+(x(3)-z1)^2);
dist2=sqrt((x(1)-y2)^2+(x(2)-x2)^2+(x(3)-z2)^2);
dist3=sqrt((x(1)-y3)^2+(x(2)-x3)^2+(x(3)-z3)^2);
sphereSurfaces contains:
function err = sphereSurfaces( x,coord,r,param )
err = [sqrt(sum((x-coord(:,1)).^2)) - r(1);
sqrt(sum((x-coord(:,2)).^2)) - r(2);
sqrt(sum((x-coord(:,3)).^2)) - r(3)];
%experimented with this instead of the above block, gives the same results, not sure which one is better/faster
% err = [sum((x-M(:,1)).^2) - r(1)^2;
% sum((x-M(:,2)).^2) - r(2)^2;
% sum((x-M(:,3)).^2) - r(3)^2];
end
The main question is above, but what about the commented block in my sphereSurfaces function? Which one should I use or is it even incorrect?
If anyone wants to run the code themselves, they can use these values. However, they were picked at pretty much random. (probably poor choice by me, test with any other values you want): y1 = 5;x1 = 5;z1 = 0;r1 = 10;y2 = 5;x2 = 5;z2 = 5;r2 = 10;y3 = 5;x3 = 7.5;z3 = 2.5;r3 = 10;

Accepted Answer

Roger Stafford
Roger Stafford on 13 Feb 2014
This is not a direct answer to your question, Mosch, but I hate to see you use the iterative methods of 'fsolve' or 'lsqnonlin' for this problem when it is possible to compute the answer directly and much more efficiently. The following is copied from the cssm newsgroup threads
http://www.mathworks.com/matlabcentral/newsreader/view_thread/239659
http://www.mathworks.com/matlabcentral/newsreader/view_thread/318366
Let p1, p2, and p3 be three-element vectors giving the x,y,z coordinates of the spheres' three centers, and r1, r2, and r3 their respective radii. Then do this:
p21 = p2-p1;
p31 = p3-p1;
c = cross(p21,p31);
c2 = sum(c.^2);
u1 = cross(((sum(p21.^2)+r1^2-r2^2)*p31 - ...
(sum(p31.^2)+r1^2-r3^2)*p21)/2,c)/c2;
v = sqrt(r1^2-sum(u1.^2))*c/sqrt(c2);
i1 = p1+u1+v;
i2 = p1+u1-v;
The i1 and i2 will be vectors for the two points of intersection. If there is any possibility of having radii such that no solution is possible (which can happen very easily) you should test that the first square root used in the calculation of v is real. An impossible problem will make it imaginary."
  2 Comments
mosch
mosch on 13 Feb 2014
Thanks! That looks like it'll solve my problem by avoiding it completely.
Giuseppe Naselli
Giuseppe Naselli on 22 Apr 2021
Hi Roger,
thanks for this.
Unfortunately the 2 url links you attached do not exist any more
However, I found your code very useful but I do not understand what each step is doing. It is any chance you can comment the lines above to explain what each one does?
Thanks in advance
G

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!