Since my days as a physics student, I have been dreaming about replacing tedious pencil and paper calculations completely by computer algebra, ideally in Python. This is the first part of a series of articles in which I will show how you can solve typical (college) math or physics problems with Python alone, specifically with sympy.

No cumbersome hand calculations are necessary! If you have an idea, how to solve the problem more elegantly with the same or other tools, please share your idea in the comments!

Today’s homework problem that we will solve with sympy is from basic analytic geometry. We have a line defined by

$$ \begin{aligned} &\vec{L}(\lambda)=\vec{a}+\lambda \vec{v}\\ \end{aligned} $$

where

$$ \begin{aligned} &\vec{a}=\left[\begin{array}{c} 3 \\ -1 \\ 0 \end{array}\right] \quad, \quad \vec{v}=\left[\begin{array}{c} 2 \\ 1 \\ -2 \end{array}\right] \end{aligned} $$

The tasks are as follows:

• Find the equation of the plane that contains the line \(\vec L\) and the point \(\vec p = (2, 1, 0)^T\)

• Find the angle between the line L and the xz-plane.

• Find the perpendicular distance between the line and the x-axis.

Finding the equation of the plane

First, start a new Jupyter notebook and define the vectors in sympy. If you don’t have sympy installed, install it by entering

pip install sympy

in the console.

For convenience, I import everything from sympy:

from sympy import *

Yes, I know this is not good style for software development because it clutters the namespace, but we are not developing software here, but solving math problems.

Now let’s define some symbols and the line:

a, lamda, mu, v, p = symbols(r"\vec{a} \lambda \mu \vec{v} \vec{p}")
L = a + lamda * v

Note that we have to write lamdawithout “d”, because lambda is a reserved keyword in Python. Oh, and yes, you can define symbols using LateX syntax like in the example above.

Now we can define the actual values of 𝑎⃗, 𝑣⃗ and 𝑝⃗ as follows:

a_ = Matrix([3, -1, 0])
v_ = Matrix([2, 1, -2])
p_ = Matrix([2, 1, 0])

Vectors are, in a certain sense, just matrices with a single column.I like to distinguish between the abstract entities 𝑎⃗ and 𝑣⃗ and the actual representations of those entities by appending a trailing underscore _. So for the line we have the representation

L_ = a_ + lamda * v_

L_

In order to proceed, we must make sure that 𝑝⃗ is not on the line, otherwise there would be infinitely many planes passing through the line and 𝑝⃗. That is easy to check, because in order for L and 𝑝⃗ to have the same z-component, we must have 𝜆=0, but then we have

L_.subs(lamda, 0)

But this is clearly not 𝑝⃗, hence 𝑝⃗ is not on 𝐿⃗. So, 𝑝⃗ and 𝑣⃗ are linearly independent and for the plane we have the parametrized form

A = a + lamda * v + mu * p

$$ \lambda \vec v + \mu \vec p + \vec a $$

I find it bit annoying in sympy (but this is similiar to other computer algebra systems like Mathematica or Maple) that it rearranges the order of the terms. I wouldn’t do this in a pencil and paper calculation, so I would love if sympy would let me keep the order. I have played around with the evaluate arguments, but that leads to other problems. So if you have an idea how to preserve order, please tell me in the comments.

Ok, we want the equation for the plane. We can get that from

$$ \left[\begin{array}{l} x \\ y \\ z \end{array}\right] \cdot \vec{n}=d $$

where 𝑛⃗ is the normal to the plane and 𝑑 is a constant. Let’s calculate the normal of the plane. This is simply the cross product of two linearly independent vectors in the plane:

n_ = v_.cross(p_)

n_

What’s the constant 𝑑? It is a constant. So choose any point in the plane and evaluate. Simply use the parametrized form for some values for 𝜆 and 𝜇. Let’s keep it simple and choose 𝜆=𝜇=0.

d_ = A.subs({a: a_, v: v_, lamda: 0, mu:0}).dot(n_)

The result is d=10. Then the equation of the plane is

from sympy.abc import x, y, z
xyz = Matrix([x, y, z])
eq = Eq(xyz.dot(n_), d_)

eq

I don’t know why, but I like to start the equation with 𝑥 and not 2𝑥. So divide the equation by 2.

eq = Eq(eq.lhs / 2, eq.rhs / 2)

eq

This is another annoyance: for manipulating equations in computer algebra systems like sympy, I have to treat left hand side and right hand side separately. I would like to write instead eq / 2 but that raises an exception. So this drives me to write a custom class inheriting from Eq:

class MyEq(Eq):
    
    def __truediv__(self, other):
        return MyEq(self.lhs / other, self.rhs / other)
    
    def __add__(self, other):
        return MyEq(self.lhs + other, self.rhs + other)


eq = MyEq(xyz.dot(n_), d_)

eq / 2

Of course, all the other operations on equations are still missing (addition, subtraction, division, power, etc.), but you get the idea.

Find the angle between the line and the (xz)-plane

xz-plane means 𝑦=0. Let’s find the angle 𝛼 between the normal to the xz-plane and the line first. The desired angle 𝛽 between the plane and the line is then 𝛽=𝜋/2−𝛼. The normal 𝑚⃗ to the xz-plane points in the y-direction, so

m_ = Matrix([0, 1, 0])

m_

For convenience, I define another function to calculate the vector norm:

def norm(vec):
    return sqrt(vec.dot(vec))

α = acos(m_.dot(v_) / norm(m_) / norm(v_))

α

So the desired angle is

β = pi/2 - α

β

Find the perpendicular distance between the line and the 𝑥-axis

Probably there are easier ways, but my first impulse is to solve this as a minimization problem. The vector joining two arbitrary points on 𝐿⃗ and the 𝑥-axis, which is parametrized by

from sympy.abc import sigma
x_axis_ = sigma * Matrix([1, 0, 0])

x_axis_

L_ - x_axis_

and the length of this vector is

dd = norm(L_ - x_axis_)

dd

Now we have to minimize this with respect to 𝜆 and 𝜎. For a minimum, the partial derivative ∂/∂𝜆 must be zero:

eq = MyEq(diff(dd, lamda), 0)

eq

I want to multiply the equation by the square root in the denominator. To do so, have a look at the arguments of the left hand side of the equation:

eq.lhs.args

(1/sqrt(4*\lambda**2 + (\lambda - 1)**2 + (2*\lambda - sigma + 3)**2),
 9*\lambda - 2*sigma + 5)

So the inverse of the square root is the zeroth argument eq.lhs.args[0]

_[0]

Now we can get rid of the ugly square root:

eq = eq / eq.lhs.args[0]

eq

Much better. Do the same with ∂/∂𝜎:

eq2 = MyEq(diff(dd, sigma), 0)

eq2

eq2 = eq2 / eq2.lhs.args[0]

eq2

The two equations for ∂/∂𝜆 and ∂/∂𝜎 must hold simultaneously, so we must solve two linear equations for the two unknowns 𝜆 and 𝜎:

result = solve([eq, eq2], (lamda, sigma))

result

{\lambda: 1/5, sigma: 17/5}

The result is a Python dict which can be used as argument in the subs method for substitution:

dd.subs(result)

And this is our result for the perpendicular distance. Let’s double check. What is the point on 𝐿⃗ where the minimum distance is obtained? It is

p_L = L_.subs(lamda, result[lamda])

p_L

And the corresponding point on the 𝑥-axis is

p_x = x_axis_.subs(sigma, result[sigma])

p_x

The two points are separated by

norm(p_L - p_x)

which agrees with our previous result.