MATLAB: Solve function not working for linear equation system with symbolic variables

solvesyms

I have set up a system of linear equations (admittedly complicated) that represent planes.
if true
  % function [S] = planefinder(N)
     P1 = [];
     P2 = [];
     P3 = [];
     T = [];
     for j = 1:size(N,1)
     P1 = [N(j,1),N(j,2),N(j,3)];
     P2 = [N(j,4),N(j,5),N(j,6)];
     P3 = [N(j,7),N(j,8),N(j,9)];
     normal = cross(P1-P2, P1-P3);
     syms x y z
     P = [x,y,z];
     planefunction = dot(normal, P-P1);
     T = [T, planefunction==0];
    end
    S = solve(T);
    end
'T' does indeed give me a system of 32 planes. However, instead of a numerical or implicit solution in terms of x,y,z I just get the following ouput:
if true
      x: [0x1 sym]
        y: [0x1 sym]
       z: [0x1 sym]
end
I am not sure what this represents; another user on a different thread with a similar problem found that a reinstall of the symbolic math package solved the problem. Is this what I need to do, and is what I've tried with solve() a method that would work? Another obvious function to try would be linsolve, but I don't know a simple way to rearrange my equations (which are of the form a*x + b*y + c*z – d ==0) so that I get a vector of the d coefficients on the right hand side. Thanks in advance.
Hi John,
I have an update on my attempt to tackle this problem. I am using a function which was written by somebody else, called ConvexPolytrope.m- attached.
This works out the convex shape that bounds half-spaces (not planes as I said before).
The input to this code as you can see is a matrix where each row is the coefficients of the half-space equation, plus an 'inner point' which the algorithm needs as a starting point.
So, in the end, I have 8 planes to which the polytrope is a solution. 2 of these define the half spaces- they are defined by my original problem.
Then there are 6 half spaces- it is simplest to choose these to be the faces of a cuboid, where one of the variables in the plane system has it's lower/upper bound.
I do this, and I am pretty sure my set of half planes is consistently defined. However, with the following input, which is supposed to be the 6 'bounding planes',
convexPolytrope([0,-1,0,-2500; 0,1,0,-2000; 0,0,-1,-50; 0,0,1,-25; -1,0,0,-2500; 1,0,0,-2000;],[0.5,0.5,0.5])
I get the same kind of output that I originally posted about
ans = [4×3 double] [4×3 double] [4×3 double] [4×3 double] [4×3 double] [4×3 double]
It's worse if I additionally supply the two planes,
convexPolytrope([1,1,-1,0; 1,0,-1,0; 0,-1,0,-2500; 0,1,0,-2000; 0,0,-1,-50; 0,0,1,-25; -1,0,0,-2500; 1,0,0,-2000;],[-0.5 -0.5 -0.5])
Error using cgprechecks (line 30) Data points containing Inf or NaN are not supported.
Error in convhulln (line 41) cgprechecks(x, nargin, cg_opt);
Error in convexPolytrope (line 38) K = convhulln(G);
What is wrong with my input that could produce the 'NaN or 'Inf' error?
    convexPolytrope.m
Matthew Worsdale

Best Answer

I'm sorry, but this is insane to do as you are doing. Solve fails because what you are doing makes no sense at all. (Do you want me to lie?) There is NO need to reinstall the symbolic toolbox.
You have sets of three points. Each set represents three points that lie in a plane (in 3 dimensions.) As it is, you then compute a normal vector using cross.
The normal vector DEFINES the plane. I'm not sure what it is that you need more than that. What are you trying to solve for?
P1 = rand(1,3)
P1 =
      0.81472      0.90579      0.12699
P2 = rand(1,3)
P2 =
      0.91338      0.63236      0.09754
P3 = rand(1,3)
P3 =
       0.2785      0.54688      0.95751
normal = cross(P1-P2, P1-P3)
normal =
     -0.23766    -0.066143     -0.18203
normal = normal/norm(normal)
normal =
      -0.7752     -0.21574     -0.59374
(There really was no need to scale normal to have unit norm, but I did. In some computations it is necessary, so this is more a matter of habit.)
The elements of normal ARE the coefficients a,b,c that form the equation of a plsne. P1 is a point in the plane (or P2 or P3.)
dot(normal,P1)
ans =
      -0.90239
So the equation of the plane that passes through those points is
syms x y z
vpa(normal*[x;y;z] == dot(normal,P1),5)
ans =
- 0.7752*x - 0.21574*y - 0.59374*z == -0.90239
(I used vpa because otherwise the coefficients are too long and hurt my eyes.)
In all of this, exactly what are you trying to solve for? There are infinitely many points in a plane. That equation defines them all.
Perhaps you are looking for a parametric form of the plane? There is NO need to use solve for that. In fact, it is a waste of CPU cycles.
nullspace = null(normal)
nullspace =
     -0.21574     -0.59374
      0.97378    -0.072158
    -0.072158      0.80142
The columns of the null space of the normal vector will span the plane.
syms u v
vpa(P1 + [u v]*nullspace',5)
ans =
[ 0.81472 - 0.59374*v - 0.21574*u, 0.97378*u - 0.072158*v + 0.90579, 0.80142*v - 0.072158*u + 0.12699]
So ANY point in the plane is described parametrically by the above expression, for some value of u and v. This is as much as solve could give you.
You comment that you are looking for the coefficients on the right hand side, thus the values of d for each plane?
As I have written it, this is no more than
d = normal*P1';
I still have no idea what it is you are trying to do. Somehow I'm not terribly confident that you do either.
Related Question