Zero division in linear equation solution

I'm trying to transform a vector to another coordinate system with different root vectors. The other root vectors are defined by three points in space that form a plane, and it's a normal vector.

{x,y,z} is the point to be transformed
{x0, y0, z0}, {x1, y1, z1} and {x2, y2, z2} are the three points
{xn, yn, zn} is the normal vector


So I solved this equation:

Solve[Apply[And,
Equal, {{x, y, z},
a ({x1, y1, z1} - {x0, y0, z0}) +
b ({x2, y2, z2} - {x0, y0, z0}) +
c {xn, yn, zn} + {x0, y0, z0}}]], {a, b, c}] // InputForm


and I got:

{{a -> -(((-(xn*y0) + xn*y2 + x0*yn - x2*yn)*(-(xn*z) + xn*z0 + x*zn -
x0*zn) - (-(xn*y) + xn*y0 + x*yn - x0*yn)*(-(xn*z0) + xn*z2 +
x0*zn - x2*zn))/((-(xn*y0) + xn*y2 + x0*yn - x2*yn)*
(-(xn*z0) + xn*z1 + x0*zn - x1*zn) -
(-(xn*y0) + xn*y1 + x0*yn - x1*yn)*(-(xn*z0) + xn*z2 + x0*zn -
x2*zn))),
b -> -((-(xn*y0*z) + xn*y1*z + x0*yn*z - x1*yn*z +
xn*y*z0 - xn*y1*z0 - x*yn*z0 + x1*yn*z0 - xn*y*z1 + xn*y0*z1 +
x*yn*z1 - x0*yn*z1 - x0*y*zn + x1*y*zn + x*y0*zn - x1*y0*zn -
x*y1*zn + x0*y1*zn)/(xn*y1*z0 - xn*y2*z0 - x1*yn*z0 + x2*yn*z0 -
xn*y0*z1 + xn*y2*z1 + x0*yn*z1 - x2*yn*z1 + xn*y0*z2 - xn*y1*z2 -
x0*yn*z2 + x1*yn*z2 + x1*y0*zn - x2*y0*zn - x0*y1*zn + x2*y1*zn +
x0*y2*zn - x1*y2*zn)),
c -> -((-(x1*y0*z) + x2*y0*z + x0*y1*z - x2*y1*z - x0*y2*z + x1*y2*z +
x1*y*z0 - x2*y*z0 - x*y1*z0 + x2*y1*z0 + x*y2*z0 - x1*y2*z0 -
x0*y*z1 + x2*y*z1 + x*y0*z1 - x2*y0*z1 - x*y2*z1 + x0*y2*z1 +
x0*y*z2 - x1*y*z2 - x*y0*z2 + x1*y0*z2 + x*y1*z2 - x0*y1*z2)/
(xn*y1*z0 - xn*y2*z0 - x1*yn*z0 + x2*yn*z0 - xn*y0*z1 + xn*y2*z1 +
x0*yn*z1 - x2*yn*z1 + xn*y0*z2 - xn*y1*z2 - x0*yn*z2 + x1*yn*z2 +
x1*y0*zn - x2*y0*zn - x0*y1*zn + x2*y1*zn + x0*y2*zn - x1*y2*zn))}}


Now if I test it with some values:

x = 0
y = 0
z = 1

x0 = -1
y0 = -1
z0 = 0

x1 = 1
y1 = -1
z1 = 0

x2 = -1
y2 = 1
z2 = 0

xn = 0
yn = 0
zn = 1

a = -(((-(xn*y0) + xn*y2 + x0*yn - x2*yn)*(-(xn*z) + xn*z0 + x*zn -
x0*zn) - (-(xn*y) + xn*y0 + x*yn - x0*yn)*(-(xn*z0) + xn*z2 +
x0*zn - x2*zn))/((-(xn*y0) + xn*y2 + x0*yn -
x2*yn)*(-(xn*z0) + xn*z1 + x0*zn - x1*zn) - (-(xn*y0) +
xn*y1 + x0*yn - x1*yn)*(-(xn*z0) + xn*z2 + x0*zn - x2*zn)))


I get:

Power::infy: Infinite expression 1/0 encountered. >>
Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>


There obviously is a solution. I mean I can calculate it in my head, it's a=0.5. Just can't figure out what I'm doing wrong.

Let's say the result of your Solve is assigned to sol. Then set the calculated value of a to be:

acalc = sol[[1, 1, 2]];


Using the particular values you suggest:

vals = {xn -> 0, x0 -> -1, y0 -> -1, z0 -> 0, x1 -> 1, y1 -> -1,
z1 -> 0, x2 -> -1, y2 -> 1, z2 -> 0, yn -> 0, zn -> 1};


does indeed lead to a problem because:

{Numerator[acalc] /. vals, Denominator[acalc] /. vals}
{0, 0}


and this is the source of your error message. If you try it for different constants (say let zn->0.1) then you get a fine solution.

One way to get around this is to set the equations to

eqns = MapThread[
Equal, {{x, y, z},
a ({x1, y1, z1} - {x0, y0, z0}) + b ({x2, y2, z2} - {x0, y0, z0}) +
c {xn, yn, zn} + {x0, y0, z0}}]


and then solve after the substitution, i.e., solve for the a, b, and c:

Solve[eqns /. vals, {a, b, c}]

{{a -> (1 + x)/2, b -> (1 + y)/2, c -> z}}


So to get x=y=z=0, you have your calculated value of 1/2.

I haven't found out why the method presented in the question doesn't produce the right formula. But if you want the formula to present an arbitrary vector that's in 3d cartesian coordinate system as a linear combination of vectors forming a basis, or in other words, want to transform the vector from 3d cartesian basis to an arbitrary one, you can do this:

The three vectors forming the basis:

v1 = {x1, y1, z1};
v2 = {x2, y2, z2};
v3 = {x3, y3, z3};


The vector in 3d cartesian coordinate system to be transformed:

p = {px, py, pz};


The math:

M = Transpose[{v1, v2, v3}];
theTransformedVector = Inverse[M].p


output:

{(pz (-x3 y2 + x2 y3))/(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 -
x2 y1 z3 + x1 y2 z3) + (
py (x3 z2 - x2 z3))/(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 -
x2 y1 z3 + x1 y2 z3) + (
px (-y3 z2 + y2 z3))/(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 -
x2 y1 z3 + x1 y2 z3), (
pz (x3 y1 - x1 y3))/(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 -
x2 y1 z3 + x1 y2 z3) + (
py (-x3 z1 + x1 z3))/(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 -
x2 y1 z3 + x1 y2 z3) + (
px (y3 z1 - y1 z3))/(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 -
x2 y1 z3 + x1 y2 z3), (
pz (-x2 y1 + x1 y2))/(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 -
x2 y1 z3 + x1 y2 z3) + (
py (x2 z1 - x1 z2))/(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 -
x2 y1 z3 + x1 y2 z3) + (
px (-y2 z1 + y1 z2))/(-x3 y2 z1 + x2 y3 z1 + x3 y1 z2 - x1 y3 z2 -
x2 y1 z3 + x1 y2 z3)}