Find plane equations given two parametric line equations

Edit: updated equation for z.

Hi,

I have two parametric line equations which intercept at (2.5,2,-2.5), I have used the below code to plot these. I believe they're perpendicular, so I am trying to work out how to find the cross product (vector normal to the two lines) and the plane equation that contains both lines.

Line (x,y,z) { x=3-t, y=1+2t, z=-4+3t } Line (x1,y1,z1) {x=9-13t, y=4-4t, z=-(5/3)-(5/3)t}

clear all;
close all;
clc;
%Parametric eq of a line
%x=x0+t*vx
%y=y0+t*vy
%z=z0+t*vy
%point
P=[3 1 -4]
v=[1 2 3];
t=-10:0.5:10;
x=P(1)-t*v(1);
y=P(2)+t*v(2);
z=P(3)+t*v(3);
x1=9-13*t;
y1=4-4*t;
z1=-(5/3)-(5/3)*t;
figure(1)
hold on
plot3(x,y,z,x1,y1,z1)
plot3(x,y,z)
hold off
grid on
xlabel('x axis')
ylabel('y axis')
zlabel('z axis')

Answer 1

[John D'Errico]

The biggest trick here is to learn to use MATLAB! MATLAB is a vector & array language. So learn to use those capabilities.

First, a line can be defined parametrically by a point on the line, and a vector that gives the direction.

Line 1:

P1 = [3 1 -4];
V1 = [-1 2 3];

So any point on line 1 is given as

P1 + t*V1

Likewise,

P2 = [9 4 5/3];
V2 = [-13 -4 -5/3];

And a point on this line is defined by

P2 + s*V2

You claim that these lines have a common point of intersection. We can test that:

Q = [2.5 2 -2.5];

So, does Q lie on line 1? If it does, then we could compute the value of t as

(Q - P1)./V1
ans =
        0.5          0.5          0.5

If Q lies on line 1, then I should expect to see THREE copies of the exact same value. So great, Q does lie on line 1. How about line 2?

(Q - P2)./V2
ans =
        0.5          0.5          2.5

So in fact, Q does NOT lie on line 2, because I do not see three copies of the same value. (The trick I used will fail only when some element of V1 or V2 is exactly zero.)

Do the lines intersect at any point? Or are they non-intersecting? Ok, I can be lazy at times. I'll let MATLAB do the work here.

syms s t
st = solve(P1 + t*V1 == P2 + s*V2,s,t)
st = 
  struct with fields:

    s: [0×1 sym]
    t: [0×1 sym]

It is not hard to show by hand either, but the symbolic toolbox makes it soooo easy. The empty result tells us that no values of s and t solve the problem. If you want to show the lines do not intersect by hand, you need to think about what an intersection would require. For example, we can view the problem in terms of linear algebra as:

A = [V1.', -V2.'];

Now, we need to solve for the 2x1 vector st, such that

A*st = P1.'-P2.'

Since the matrix A is a 3x2 matrix, no exact solution will exist, unless this rank is 2:

rank([A,P1.'-P2.'])
ans =
     3

So the lines are indeed skew and non-intersecting. Were V1 and V2 the same vector (or a simple multiple of each other) then the lines would be parallel. In this case of course, we again see no intersection exists.

So the very first premise of your question appears to be false, that the lines do intersect. But if the lines are indeed non-parallel and non-intersecting, then you cannot compute a plane that contains them.

Oh, one other point, V1 and V2 are indeed perpendicular vectors, so you got that right.

dot(V1,V2)
ans =
     0

My guess is your biggest problem is the definition of your two lines. I think you have z written down wrong for line 2.

Finally, at the point where you will need to determine the plane equations, there are several tools that will work. I would suggest either null or cross as good tools here, rather than writing the cross product computations explicitly. But first, you need to resolve the issue I brought up.

[FortuitousMonkey]

Thank you for the really informative feedback. You're correct I did make a mistake with Z the '-' was missing from V2(1) and actually created V2 and P2 I've smartened up my code now as well. Though I didn't use matlab to arrange or solve the equations I did that by hand (still learning as you can see). You have also helped me understand deriving the points.

if true
P2 = [9 4 -5/3];
V2 = [-13 -4 -5/3];
%Equations
x1=P2(1)+t*V2(1);
y1=P2(2)+t*V2(2);
z1=P2(3)+t*V2(3);

I checked the dot of V1 V2 and found 0, so they're still perpendicular.

dot(V1,V2)
ans =
   0

I then checked if Q lies on both lines it does. From my understanding the output is t when the intercept =Q, correct?

Line1

    (Q - P1)./V1
ans =

    0.5000    0.5000    0.5000

Line 2

(Q - P2)./V2
ans =
  0.5000    0.5000    0.5000

I then ran your final section of code, and found the output had changed. Though I will say I don't understand fully what the output st is describing.

syms s t
st = solve(P1 + t*V1 == P2 + s*V2,s,t)
st = 
struct with fields:

    s: [1×1 sym]
    t: [1×1 sym]

If my understanding is correct I now have 2 perpendicular lines which intersect at Q and we can proceed to look at deriving the plane equation?

I've added a plot should it be helpful for someone elses understanding.

[John D'Errico]

Good. I had a funny feeling the problem was in z for line 2.

st = solve(P1 + t*V1 == P2 + s*V2,s,t)
st = 
  struct with fields:

      s: [1×1 sym]
      t: [1×1 sym]

st.s
ans =
1/2

st.t
ans =
1/2

so the solution happens when s=t=1/2.

Are the lines orthogonal? I.e., perpendicular? dot tells us they are so, but they don't look as if they are. This is a graphics issue that befalls many of us. It is easy to trip over.

T = linspace(-5,5,100)';
S = linspace(-5,5,100)';

L1 = P1 + S.*V1;
L2 = P2 + T.*V2;

plot3(L1(:,1),L1(:,2),L1(:,3),'r-')
hold on
plot3(L2(:,1),L2(:,2),L2(:,3),'b-')
box on
grid on

Hmm. No matter how I rotate that plot around, I cannot get the lines to LOOK as if they are perpendicular.

The trick is to make the axes have the same effective units.

axis equal

Now we can see they are indeed perpendicular.

Ok. Now, what is the definition of a plane that contains these two lines? How will we write it?

For any point X where X is a vector that lies in R^3, so the 3-dimensional space where these lines live,

dot(Q - X,N) = 0

Here Q is a point in the plane. I've chosen the point of intersection, but any point on that plane would suffice. N is the normal vector to the plane.

N = cross(V1,V2)
N =
       8.6667      -40.667           30

N = N/norm(N)
N =
      0.16903     -0.79314      0.58511

So cross produces the vector orthogonal to both V1 and V2, though the result is not a unit vector. Dividing by the norm solves that. I have done so here mainly to show that cross and null would produce the same result (except for possibly a change in sign.)

N = null([V1;V2])
N =
      0.16903
     -0.79314
      0.58511

The equation of your plane is given by the dot product above.

syms x y z
(Q - [x,y,z])*N == 0
ans =
(61*35^(1/2)*(y - 2))/455 - (35^(1/2)*(x - 5/2))/35 - (9*5^(1/2)*7^(1/2)*(z + 5/2))/91 == 0

In terms of decimal numbers, the plane equation is:

vpa(ans,10)
ans =
0.7931447621*y - 0.1690308509*x - 0.5851067917*z - 2.626479376 == 0.0

We can go back and now confirm that both lines fall in this plane.

syms s
vpa((Q - (P1 + s*V1))*N)
ans =
0

vpa((Q - (P2 + s*V2))*N)
ans =
2.7550648847397363468017262591146e-40 - 5.5101297694794726936034525182293e-40*s

We nailed it for the first line. I hope you will forgive MATLAB for the tiny error down in the 40th decimal place for line 2. ;-)