Mathcad Solution To XKCD Problem

Quote of the Day

I don't need time – I need a deadline.

— Duke Ellington, describing how he is motivated to get something done. I must admit that I need deadlines for writing, otherwise I would never get anything written.

Introduction

Figure 1: Velociraptor Vs Prey Problem Statement (Source).

XKCD is a great comic strip by Randall Munroe that takes a quirky look at the world of science. Randall recently posted a set of questions for a substitute teacher to pose to a math class that were interesting and attracted the attention of some problem-solving enthusiasts.

I thought I would jump into the fray by solving the problem using Mathcad and also taking my first dive into Python.

Background

While this is a simple problem, it is interesting to see how different people approach it. Here are two other solutions that I recommend you read:

• Jon Peltier's Excel Solution

  Jon is a world-class Excel expert (e.g. MVP) with a real gift for Excel charting. He attacks the problem using a standard spreadsheet approach and a Visual Basic for Applications (VBA) approach.

• Rhett Allain's Python Solution

  Rhett wrote a Python-based solution that I used as the basis for my first Python program. I modified his routine to use trapezoidal integration instead of Euler, and I implemented a different scheme for the clipping the velocities.

Analysis

Mathcad Solution

Figure 2 shows how I used Mathcad to solve the substitute teacher problem.

Figure 2: Mathcad Solution to XKCD Problem.

Figure 3 shows my graphical view of the solution.

Figure 3: Plot of Raptor Versus Prey.

Python Solution

I have never written a Python program before and this seemed like a great opportunity to learn. When I looked at Rhett's program, it reminded me a simplified form of Java – my usual programming language choice. So I quick ran through the Python course on Code Academy and loaded Python onto my Eclipse install.

Figure 4 shows my Python routine, which gave me the same results as Mathcad. Python seems a lot less wordy than Java – what were the Java folks thinking?

import matplotlib.pyplot as plt
#!/usr/bin/env python
"""
Abstract: 
This module iteratively solves the xkcd substitute teacher problem:

    The velociraptor spots you 40 meters away and attacks, accelerating
    at 4 m/s^2 up to its top speed of 25 m/s. When it spots you, you 
    begin to flee, quickly reaching your top speed of 6 m/s. How far can
    you get before you are devoured.

My intent here is to solve the problem with a simple Python program. My 
approach here is a slightly modified version of Rhett Allain's Python code 
at http://www.wired.com/2015/10/heres-solve-xkcd-velociraptor-problem-code/
My modifications are simple. 
1. I used trapezoidal rule rather than Euler method for the ODE solution. 
2. I decided to store the intermediate results in lists.
3. Implemented the clipped velocities using minimum function.

Rhett assumes that a prey will accelerate at 3 m/s^2 on their way to a top
speed of 6 m/s, so I will too.
"""
# Initial Conditions
xv=-40      #this is the initial location of velociraptor
xp=0        #location of the prey
av=4        #acceleration of velociraptor
ap=3        #accel of prey 
vvmax=25    #maximum velocity of the velociraptor
vpmax=6     #max velocity of prey
vp=0        #starting velocity of prey
vv=0        #starting velocity of velociraptor
t=0         #starting time
dt=0.0001   #time step
tlist=[]                #Store all time values computed
xpreylist=[]            #Store all prey positions
xvelociraptorlist=[]    #Store all raptor positions

#this is the loop. It runs until the velociraptor
#catches up to the prey
while xv<=xp:#first check if the prey is at max v
   
    #calc new prey velocity after time interval
    vp0=vp
    vp=min(vp+ap*dt,vpmax)
  
    #calc new velociraptor velocity
    vv0=vv
    vv=min(vv+av*dt,vvmax)
    
    #calc new positions
    xp=xp+(vp+vp0)*0.5*dt   #Trapezoidal integration
    xpreylist.insert(0,xp)
    xv=xv+(vv+vv0)*0.5*dt   #Trapezoidal integration
    xvelociraptorlist.insert(0,xv)
    tlist.insert(0,t)
    #update time
    t=t+dt

print 'Prey Position @ intercept time: {0:4}'.format(xp)
print 'Intercept Time (s): {0:4}'.format(t)

# Plot the prey and velociraptor positions

plt.plot(tlist,xpreylist)
plt.plot(tlist,xvelociraptorlist)
plt.ylabel('Position (m)')
plt.xlabel('Time (s)')
plt.suptitle('Velociraptor Vs Prey Position', fontsize=20)
plt.show()

Figure 4: My Python Solution for the Substitute Teacher Problem.

Python Graphic Output

Figure 5 shows the output from my Python program. The results are the same as given by Mathcad.

Figure 5: Python Plot of Prey and Velociraptor Position.

Conclusion

This problem nicely illustrates how different tools can be used to solve a problem. In my case, I tend to use computer-algebra systems, like Mathcad and Mathematica, to experiment with different approaches to solve a problem. I only commit to software when I need more speed or the ability to handle large scale.