Simple calculations giving high imprecisions

I'm doing the following simple code, that aims to find a point position knowing distances in 2D

but the outputs is far from the actual values:

the gap for x is around 0.4! which is huge, why is it so? normally the code is correct and one solution should match almost exactly the actual position

x =

-22.5971  -70.9058

x3 =

-22.1513

y =

   34.5832   35.4659

y3 =

   34.5751

----code---

X =[ -47.0124   47.5995;
  -47.6018   15.3411;
  -22.1513   34.5751];

x1 = X(1,1);
y1 = X(1,2);
x2 = X(2,1);
y2 = X(2,2);
x3 = X(3,1);
y3 = X(3,2);

r1 = norm(X(3,:) - X(1,:));
r2 = norm(X(3,:) - X(2,:));

B = (sum(X(2,:).^2)-sum(X(1,:).^2)-r2^2+r1^2)/2;

% trying to retrieve X(3,:) knowing the distances only

% (x3-x1)^2 + (y3-y1)^2 =r1^2  (1)
% (x3-x2)^2 + (y3-y2)^2 =r2^2  (2)
% (2) - (1) to have a relation between x3 and y3
% x3(x2-x1) + y3(y2-y1) = B
% we replace in (1)

if y2-y1~=0
  'we should be there'

    C = y1-B/(y2-y1);
    tau = (x2-x1)/(y2-y1);

    a = 1+tau^2;
    b = -2*(x1+ tau*C);
    c = x1^2+C^2-r1^2;

    delta = b^2-4*a*c

    x = (sqrt(delta)*[1 -1] - b)/(2*a) % the 2 solutions for x

    x3 % check actual value

    y = (B-x*(x2-x1)) / (y2-y1)

    y3 % actual value

end

Answer 1

[Geoff Hayes]

cyril - I think that there is a sign mismatch in the calculation of C. Rather than

 C = y1-B/(y2-y1);

it should be

 C = -y1+B/(y2-y1);

due to the re-arrangement of

 x3(x2-x1) + y3(y2-y1) = B
 => y3 = [B - x3(x2-x1)]/(y2-y1)
 => y3 = B/(y2-y1) - x3(x2-x1)/(y2-y1)

and the above substitution into

 (x3-x1)^2 + (y3-y1)^2 =r1^2 
 => (x3-x1)^2 + (B/(y2-y1) - x3(x2-x1)/(y2-y1) - y1)^2 = r1^2
 => (x3-x1)^2 + (C - x3*tau)^2

where

 C   = B/(y2-y1) - y1
 tau = (x2-x1)/(y2-y1)

Note that a couple of other things you could do is to replace the norm calls with just

 r1 = (x3-x1)^2 + (y3-y1)^2; %norm(X(3,:) - X(1,:));
 r2 = (x3-x2)^2 + (y3-y2)^2; %norm(X(3,:) - X(2,:));

to avoid the square root of norm and subsequent squaring (where needed) of r1 and r2. (Note that if you use the above two replacements, then you will have to remove the squares for r1 and r2.)

When I re-run your code with the above changes, I get the expected answer:

 x =
  -22.1513  -72.3327

 x3 =
  -22.1513

 y =
   34.5751   35.4920

 y3 =
   34.5751

[cyril]

thanks very very much Geoff for the error