# The XKCD Raptor Problem¶

This is originally from an XKCD comic https://xkcd.com/135/

The problem is to figure out if a human at the center of an equilateral triangle 20 m on each side will survive if a raptor is placed at each of the vertices. One of the raptors is wounded.

In this sheet, we simulate a mathematical model of the problem that we built up in class.

• The human runs at 6.7 m/s in a straight line.
• A raptor runs at ~40 mph, or 17 m/s.
• A wounded raptor runs at about 10 m/s.
• Raptors run directly at the humans at their maximum velocity.
• We neglect the impact of acceleration.
• We assume that raptors can catch a human within a few centimeters.
• We assume there are no limits to the raptor motion in chasing the human.
• We assume the human only runs in a single angle.

The equations of motion for a human are then simple: $$\frac{dh(t)}{dt} = \begin{bmatrix} \cos \theta \\ \sin \theta \end{bmatrix} v_{\text{human}}$$

The raptors chase the human $$\frac{dr(t)}{dt} = \frac{h(t) - r(t)}{\| h(t) - r(t) \|} v_{\text{raptor}}$$

Note that we will have three positions for the raptors, one for each. There will also be a special velocity for the wounded raptor.

To make the problem slightly easier, we position the raptors 20m away from the human.

In [2]:
angle=0.
vhuman=6.7
vraptor0=10. # the slow raptor velocity in m/s
vraptor=17. #

raptor_distance = 20.

raptor_min_distance = 0.2; # a raptor within 20 cm can attack
tmax=10. # the maximum time in seconds
nsteps=1000

"""
This function will compute the derivatives of the
positions of the human and the raptors
"""
function compute_derivatives(h,r0,r1,r2)
dh = [cos(angle),sin(angle)]*vhuman
dr0 = (h-r0)/norm(h-r0)*vraptor0
dr1 = (h-r1)/norm(h-r1)*vraptor
dr2 = (h-r2)/norm(h-r2)*vraptor
return dh, dr0, dr1, dr2
end
"""
This function will use forward Euler to simulate the Raptors
"""
function simulate_raptors()
# initial positions
h = [0.0,0.0]
r0 = [1.0,0.0]*raptor_distance
r1 = [-0.5,sqrt(3.)/2.]*raptor_distance
r2 = [-0.5,-sqrt(3.)/2.]*raptor_distance

# how much time el
dt = tmax/nsteps
t = 0.

hhist = zeros(2,nsteps+1)
r0hist = zeros(2,nsteps+1)
r1hist = zeros(2,nsteps+2)
r2hist = zeros(2,nsteps+2)

hhist[:,1] = h
r0hist[:,1] = r0
r1hist[:,1] = r1
r2hist[:,1] = r2

for i=1:nsteps
dh, dr0, dr1, dr2 = compute_derivatives(h,r0,r1,r2)
h += dh*dt
r0 += dr0*dt
r1 += dr1*dt
r2 += dr2*dt
t += dt

hhist[:,i+1] = h
r0hist[:,i+1] = r0
r1hist[:,i+1] = r1
r2hist[:,i+1] = r2

if norm(r0-h) <= raptor_min_distance ||
norm(r1-h) <= raptor_min_distance ||
norm(r2-h) <= raptor_min_distance
@printf("The raptors caught the human in %f seconds\n", t)

# truncate the history
hhist = hhist[:,1:i+1]
r0hist = r0hist[:,1:i+1]
r1hist = r1hist[:,1:i+1]
r2hist = r2hist[:,1:i+1]

break
end
end
return hhist, r0hist, r1hist, r2hist
end

angle = pi/4.
simulate_raptors();
The raptors caught the human in 1.240000 seconds
In [3]:
using Plots
In [4]:
function show_raptors()
hhist, r0hist, r1hist, r2hist = simulate_raptors()
plot(vec(hhist[1,:]),vec(hhist[2,:]),linewidth=3)
plot!(vec(r0hist[1,:]),vec(r0hist[2,:]),color=:red)
plot!(vec(r1hist[1,:]),vec(r1hist[2,:]),color=:red)
plot!(vec(r2hist[1,:]),vec(r2hist[2,:]),color=:red)
plot!(grid=false,xticks=false,yticks=false,legend=false,border=false)
plot!(xlim=[-20.,20.],ylim=[-20.,20.])
annotate!(r0hist[1,1], r0hist[2,1]-1,text("Wounded Raptor",:right,"cmss10"))
annotate!(r1hist[1,1], r1hist[2,1],text("Raptor 1",:right, "cmss10"))
annotate!(r2hist[1,1], r2hist[2,1],text("Raptor 2",:right, "cmss10"))
annotate!(hhist[1,end], hhist[2,end],text(" Crunch",:left,:red, "cmss10"))
end
angle = pi/4
show_raptors()
The raptors caught the human in 1.240000 seconds
[Plots.jl] Initializing backend: pyplot
WARNING: No working GUI backend found for matplotlib.
Out[4]:
In [5]:
angle = pi/7.
show_raptors()
The raptors caught the human in 1.390000 seconds
Out[5]:
In [35]:
nsteps=10000
angle = pi/7.
show_raptors()
The raptors caught the human in 1.381000 seconds
Out[35]:
In [41]:
raptor_min_distance=0.00001
println("This is non-physical")
show_raptors()
This is non-physical
Out[41]: