python calculus

9 min readJan 8, 2021

(part 2) “the goods”

SymPy works, let’s do something with it.

I’ve opened my old calculus text book and have found a classic calc 3 problem, the classic “three dimensional motion” problem, where you can infer a parametric equation for acceleration, and you integrate on through to get to a position function! check it out:

The problem wants you to find the distance the projectile has travelled in meters.

OKAY… our goal:

we want to write a python function that will take in the parameters given (the angle above the horizontal in degrees, the initial velocity of the projectile right after being fired in meters per second and the y component of acceleration in meters per second squared) and will return the distance the projectile has gone in meters.

First

I am going to do the problem on paper so that, A. we understand the problem and have a firm grasp of whats going on at each step, and B. we can check the inputs/outputs at each step and can thus ensure our function is behaving the way we want it to.

…deep breaths, and:

what we’ve done:

1.We have enough data to assume an acceleration function, call it a(t). We know the projectile is accelerating down towards the earth at a rate of 9.8 meters per second squared (this will be expressed as a negative z component as down is negative z). And, we know that the projectile accelerates north at a rate of 0.36 meters per second squared (positive y, which in our model is north), so:

a(t) = < 0, 0.36, -9.8 >

2. We can then integrate this to get the velocity function (the derivative of velocity is acceleration, so let’s go backwards, integrate!) Oh, and we want to do this parametrically! (which just means that each component will change with respect to ONE variable, we will make that variable t (for time). So we integrate each component with respect to t:

∫a(t) = v(t) = < 0, 0.36t, -9.8t > + 𝘾

3. This gives us the (the velocity function + 𝘾); 𝘾 being the arbitrary constant. We need to find 𝘾 in order to find the true velocity curve, so we look at what is happening at the moment it is fired (when t = 0, for the velocity curve). The third piece of data we were given was the at time zero the projectile was going 300 meters per second ( |𝙑₀|=300 m/s ). Remember, with derivative vector functions, this means that the magnitude of this vector is 300 (magnitude = speed). So we can turn this into a little triangle and use pythag! as in the little picture I drew:

This represents the vector produced by the vector valued function of the velocity curve at time zero. We can now find the arbitrary constant (the 𝘾 being added to the x and z components, as this little triangle is situated between the x and z axis). So to describe the triangle above:

<300cos(30), 0, 300sin(30)>
simplifies to => <150√(30), 0, 150>, this is 𝘾!

we add this 𝘾 to the vector valued function we got earlier in order to get the true velocity function:

<150√(30), 0, 150> + < 0, 0.36t, -9.8t >Our Velocity function => v(t) = <150√(30), 0.36t, -9.8t + 150 >

3. We want a position function r(t), and now that we have Velocity, we can do that (the derivative of position is velocity is position). So we do it again, integrate:

∫v(t) = r(t) = <150√(30)t, 0.18t², -4.9t² + 150t> + 𝘾

Since our model assumes that we fired the projectile from the origin, we can just assume that our 𝘾 is the “zero vector” or: <0,0,0>. Thus, our position function is:

Position function => r(t) = <150√(30)t, 0.18t², -4.9t² + 150t>

4. Now that we understand the position of our projectile as it moves through time, we can find when it hits the ground. Our goal is to “find the distance it travelled”, and in order to do so we need to know when it left the ground and when it hit the ground. From there, we can find where it left the ground and where it hit the ground, and get the distance. Since “up” and “down” are represented by our Z axis, the “ground” is zero on the Z axis…. thus if we find where our Z component equals zero in terms of time, we will know when the projectile is in the air and when it’s on the ground. So we solve our Z component for the value of t that makes it equal to 0:

-4.9t² + 150t = 0
4.9t² = 150t
4.9t = 150
t = 150/4.9 <= this, in seconds, is when our projectile hits the    ground!

(note* this is a quadratic… because the “ t² ”, so that means there are 2 values where the equation equals 0. If we weren’t 100% sure that this curve starts at the origin and ends somewhere else we would have to look at the values and throw out non-sensical ones. Like, say the curve equals 0 at a negative t value…we can’t have “negative time” so we throw that one out).

Now that we know how many seconds after the projectile was fired it took for it to hit the ground, we can plug that value of time into our X and Y components of our position function in order to find out “how much X and Y” the projectile moved:

t = 150/4.9r(t) = <150√(30)t, 0.18t², -4.9t² + 150t>r(150/4.9) = <150√(30)(150/4.9), 0.18(150/4.9)²>solving both components:Xdistance ≈ 7953.294525
Ydistance ≈ 168.6797168

5. We have the amount our projectile has moved in the X and Y directions, see my little picture:

Just another triangle! So to find the total distance, we can just use o’l faithful, pythag:

totalDistance = √(7953.294525)² + (168.6797168)²totalDistance ≈ 7955.08307m

The total distance the projectile has travelled is about 7955.08307 meters…

NOW LETS WRITE IT IN PYTHON:

Now that we understand the problem, let’s begin!

In calc.py, lets:

import SymPy
define the variables we are going to be working with as Symbols, we only need time (t)
define the our given parameters
make our function (call it find_distance ), have it print something to make sure it is working

Okay cool, we’re getting a tuple back containing our 3 pieces of data… First order of business SymPy likes radians. The angle of launch we were given is in degrees, so let’s convert:

180° = 𝜋
=> 1°= 𝜋/180
We were given 30° so:
30°𝜋/180
simplifies to => 𝜋/6

If our function is one that will take in degrees, multiply them by 𝜋/180 and we have radians. Lets test it (we should get back 𝜋/6). Also, this function is going to be working with the constant for acceleration of gravity at the earth’s surface, so lets define that constant as well, call it gC (gravitational constant):

Yay, we converted to radians.

Since we are going to be working with Vectors and Vector valued functions, we need to find a good way to do this. We are going to use SymPy’s CoordSys3d, which allows you to import vectors as a class and different manipulations of those vectors are instances of different “vector subclasses” (you can read about it here and here). We will be able to define a 3d Cartesian coordinate system, and make instances of i j k style vectors (you can read about those here). We start by making and naming a coordinate system, and then defining a vector reaching from the origin to that (i , j, k) point. If the vectors value along a particular axis is 0, we just don’t write that component. For instance if we have a vector that starts at the origin and reaches to the point (1,3,0), we write that as: i +3k:

from sympy.vector import CoordSys3D
from sympy import *N = CoordSys3D('N')myVector = N.i + 10*N.kmyVector is a Vector pointing from the origin to the point (1,0,10)
(N.i means i)

We can also multiply in variables, and can then make parametric equations! So, we can make our acceleration function:

Now we can integrate to get our v(t) + 𝘾… Which, in SymPy is as easy as it sounds! All you do is:

myFunction = integrate(Function, x) 
"integrate the function called "function" with respect to x"

let’s integrate with respect to t! (note* we can add in SymPy’s .evalf() to get decimal approximations as well as its .doit() method which “Evaluate objects that are not evaluated by default like limits, integrals, sums and products”… remember, everything is “symbolic” and is kept exact unless you tell SymPy otherwise. This can help with the process of checking values (you can read about .doit() and .evalf() here)

and lets check:

recall that our velocity function before we found 𝘾 was:
v(t) = < 0, 0.36t, -9.8t > + 𝘾

So we are good so far!

Lets get that 𝘾:

Our 𝘾 was <150√(30), 0, 150>, lets check:

Now for our true velocity function:

before it was: v(t) = <150√(30), 0.36t, -9.8t + 150 >… lets check:

Awesome…. Now we integrate for position (remember that this 𝘾 was 0).

Lets compare to our position function which was:
r(t) = <150√(30)t, 0.18t², -4.9t² + 150t>

cool.

Now we need to grab this Z component and solve for when t = 0, we’ll call this the “zeroTimes”. Pictured below, on lines 5, 6 and 7 we are going into our position function, extracting the -4.9t² + 150t, and making it so that we can work with it as a solvable expression. On line 4 we solve that expression for the values that make it equal to 0:

Remember because this is t² there are extraneous solutions… SymPy will return us 2 values, see:

Our seconds are correct! (30.61224489… and 150/4.9 are the same) Only problem is that we only want the landing time, so let’s filter the tuple and grab the only positive time. (Again I’m taking liberties with the edge cases, I’m going very very specific, and we can later refine this to accept more and more cases… but for now).

Now we grab that value and throw it into the X and Y components in order to get the X and Y distances:

above we, grab the X and Y components from the position function, and solve both for the landing time (.subs() means substitute). Our original distances were:

Xdistance ≈ 7953.294525
Ydistance ≈ 168.6797168