Hi John,
What the solver is trying to tell you is that a solution does not exist for this problem, or more correctly it exists only in limited circumstances. Suppose the three values contained in the equation expressions are denoted by a column vector q. In your case,
q = [-2.9404; 2.1632; 2.3319]
For an arbitrary q there is no solution for p. If q is restricted in the manner mentioned below, then there are an unlimited number of solutions for p.
The three equations show that after the matrix multiplication, you are interested in the contents of the last column. For any 4x4 matrix G the last column is
G*[0;0;0;1] % [0;0;0;1] is a column vector
In this case it is
(A*B*C)*[0;0;0;1]
where A,B,C are your three matrices in easier notation. To simplify things you can multiply C*[0;0;0;1] to obtain a column vector, then multiply that result by B to obtain another column vector, etc. The net result is the column vector
[(-R*p +p);0]
and the desired solution is
(-R*p +p) = q
(I-R)*p = q % I = 3x3 identity matrix
p = inv(I-R)*q % this woud be (I-R)\q in practice
However R, being a rotation matrix, has one eigenvalue equal to 1. Subtracting constant*I off of a matrix reduces all its eigenvalues by the same amount, meaning that (I-R) has an eigenvalue of 0. The inverse of (I-R) does not exist. Hence the huge values when trying to calculate p.
There is a restriction on q that does allow a solution. Suppose u is the eigenvector corresponding to eigenvalue 0 of (I-R). If q has no component that is parallel to u, then a set of solutions does exist. From the given R and q there are some indications that that was the original intent.
Your rotation matrix must actually have a lot of significant figures but is only shown to 5 significant figures. The smallest eigenvalue of (I-R) is about -3.68e-5, but to illustrate the result, it's necessary to have a more accurate R. The following code creates an R very close to the first R, but with smallest eigenvalue -6.2427e-17.
R5 = [0.8769 -0.0366 0.4793
0.2583 0.8769 -0.4055
- 0.4055 0.4793 0.7784];
D = R5(:,1);
E = R5(:,2);
E = E/norm(E);
D = D-(D'*E)*E; % make D perpendicular to E (Gram-Schmidt)
D = D/norm(D);
F = null([D E]');
R = [D E F]
If [v lam] = eig(I-R), then v consists of the eigenvectors in columns, with v(:,3), as it turns out, corresponding to eigenvalue 0. In what follows let
q = v*qbar qbar = v\q p = v*pbar pbar = v\p
with q, qbar, p, pbar being column vectors. So qbar consists of the amplitudes of the three eigenvectors making up q, and the required condition on q is qbar(3) = 0:
I = eye(3,3);
[v e] = eig((I-R),'vector')
qbar = v\q
qbar(3)
ans = 1.7028e-04 - 1.5283e-16i % not wanted
qbar(3) = 0
qnew = v*qbar;
Now it's possible to invert (I-R) in the reduced vector space:
p = (I-R)\qnew
p =
-7.9119
-10.0159
3.3371
(I-R)*p -qnew % should be small (zero by construction)
So that's a solution. However, (I-R)*v(:,3) = 0, meaning that (I-R) kills off any contribution of v(:,3) within p. That means that if p is a solution, then
p + const*v(:,3)
is also a solution. Therefore pbar(3) can be an arbitrary constant:
pbar = v\p
pbar(3) = 100; % for example
pnew = v*pbar
pnew =
68.8738
66.7815
28.9263
(I-R)*pnew-qnew % should be small
ans =
1.0e-14 *
-0.5773
-0.2220
0.3997
As you can see the new p is quite different from the original one. In practice if pbar(3) is as large as +-1e6, the solution is still good to a part in 10^10.
