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 m/s in a straight line.
A raptor runs at ~40 mph, or 15 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 [22]:

    
vhuman=6.0
vraptor0=10.0 # the slow raptor velocity in m/s
vraptor=15.0 # 

raptor_distance = 20.0

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(angle,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(angle)
    # 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(angle,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(angle);









    



WARNING: Method definition compute_derivatives(Any, Any, Any, Any, Any) in module Main at In[20]:16 overwritten at In[22]:16.
WARNING: replacing docs for 'compute_derivatives :: Tuple{Any,Any,Any,Any,Any}' in module 'Main'.
WARNING: Method definition simulate_raptors(Any) in module Main at In[20]:27 overwritten at In[22]:27.
WARNING: replacing docs for 'simulate_raptors :: Tuple{Any}' in module 'Main'.






    



The raptors caught the human in 1.410000 seconds



In [3]:

    
using Plots



In [29]:

    
function show_raptors(angle)
    hhist, r0hist, r1hist, r2hist = simulate_raptors(angle)
    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))
    annotate!(r1hist[1,1], r1hist[2,1],text("Raptor 1",:right))
    annotate!(r2hist[1,1], r2hist[2,1],text("Raptor 2",:right))    
    annotate!(hhist[1,end], hhist[2,end],text(" Crunch",16,:red,:left))    
    title!(@sprintf("Survival time = %f sec",(-1+length(hhist[2,:]))*tmax/nsteps))
end
angle = pi/4
show_raptors(angle)









    



The raptors caught the human in 1.410000 seconds






    



WARNING: Method definition show_raptors(Any) in module Main at In[28]:2 overwritten at In[29]:2.






    Out[29]:



In [30]:

    
angle = pi/7.
show_raptors(angle)









    



The raptors caught the human in 1.400000 seconds






    Out[30]:



In [18]:

    
nsteps=10000
show_raptors(pi/7.)









    



The raptors caught the human in 0.854000 seconds






    Out[18]:



In [31]:

    
using Interact
@manipulate for angle=0:0.01:2*pi
    show_raptors(angle)
end









    











    



The raptors caught the human in 1.260000 seconds






    Out[31]:



In [ ]:



In [19]:

    
raptor_min_distance=0.00001
println("This is non-physical")
show_raptors(angle)









    



This is non-physical






    Out[19]: