Solve function not working for linear equation system with symbolic variables

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 = [4x3 double] [4x3 double] [4x3 double] [4x3 double] [4x3 double] [4x3 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

Answer 1

[John D'Errico]

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.

Answer 2

[Matthew Worsdale]

Hello John

Thanks for your reply. Sorry I should have explained my motivation and I suspect I am working under some false assumptions here. Essentially, my task is to find the convex polytrope which bounds a set of planes in 3D. MATLAB has in the past, I think, supplied a function for this, which I am attaching. (but it may have not been published by MATLAB, I am not sure). This function takes a set of plane equations as an input, plus an inner point where a number of the planes cross.

The inner point is what I am trying to obtain (I have the equations of the planes) I started out with 8 points in 3D, where my 'z' values are determined as a function of x,y.

I have minimum and maximum x,y values, so there are 4 pairs of coordinates and 2 z values that each corresponds to.

I worked out the equations of all the planes you could make passing through 3 of the points. This is why I use planefunction- my function 'planefinder' takes a list of point triplets I chose to make the planes.

That gives 32 planes in total! You could work with a smaller number, and I suspect I will have to given the complexity, but didn't know how to choose the relevant subset that is needed for the solution (if indeed it exists). So to summarise- I am looking for intersection point(s) of my set of planes.

My approach was any function of MATLAB that can solve a large set of simultaneous equations- I guess I need to produce a smaller set.

Matt

[John D'Errico]

Ok. so you have 32 planes that you want to solve for. Are you sure it is only 32? nchoosek(8,3) = 56. So there are 56 possible ways to choose 3 out of 8 points, thus 56 possible planes given 8 points.

Regardless, I'm not sure what you mean by a "convex polytope that bounds a set of planes". Planes are infinite in extent.

Part of me is wondering if your goal is simply the convex hull of your original set of 8 points. So please explain more. I'm pretty sure there is a simple solution to this, once I understand your goal.

[Matthew Worsdale]

Hi John,

I did not see this comment.

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 = [4x3 double] [4x3 double] [4x3 double] [4x3 double] [4x3 double] [4x3 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?