Julia Community 🟣

Gage Bonner
Gage Bonner

Posted on

Firecrackers and Envelopes


At the risk of going full recipe blog, I'm going to start with a boring but harmless personal anecdote.

When I was a wee lad, I thought that the secret to fame and fortune was solving Project Euler problems. If you aren't familiar, they describe it as,

Project Euler is a series of challenging mathematical/computer programming problems that will require more than just mathematical insights to solve. Although mathematics will help you arrive at elegant and efficient methods, the use of a computer and programming skills will be required to solve most problems.

It would still be a couple years before I discovered Julia, but armed with my rudimentary understanding of CALC101 and even more rudimentary understanding of Python 2.7, I got to work. I'm neither famous nor rich so the obvious conclusion is that I didn't solve enough problems, but I hacked my way through a few. Eventually it became clear that I simply didn't have the algorithmic chops (or patience) to make much more headway so I transitioned from solving problems to just sort of skimming them and seeing if I could think of anything in five minutes. You might be surprised to learn that this actually worked for one (1) problem.

I stumbled onto Problem 317: Firecracker. There were no gigantic prime numbers or arrays of obtusely defined integers, it was just good old projectile motion.

A firecracker explodes at a height of 100m100\,\text{m} above level ground. It breaks into a large number of very small fragments, which move in every direction; all of them have the same initial velocity of 20m/s20\,\text{m/s} .

We assume that the fragments move without air resistance, in a uniform gravitational field with 9.81m/s29.81\,\text{m/s}^2 .

Find the volume (in m3\text{m}^3 ) of the region through which the fragments move before reaching the ground. Give your answer rounded to four decimal places.

This problem will be today's topic. I'll show you the solution I came up with all those years ago, and atone for my past sins by also providing a Julia solution. For complete transparency, reading further will spoil the problem!

The Codeless Solution

First, clearly we need some physics - without knowing how a projectile behaves this problem is simply impossible (e.g. what if it had air drag?) A quick skim of the first chapter of any physics textbook will tell us the height of the firecracker yy as a function of its horizontal displacement xx is given by

y=y0+xtanθgsec2θ2v2x2 y = y_0 + x \tan \theta - \frac{g \sec^2 \theta}{2 v^2} x^2

where y0=100my_0 = 100\,\text{m} is the initial height, v=20m/sv = 20\,\text{m/s} is the initial velocity, g=9.81m/s2g = 9.81\,\text{m/s}^2 is the acceleration due to gravity and θ\theta is the launch angle with respect to the horizontal. Let's get plotting; we'll assume that the firecracker is launched to the "right" and only look at the part of the trajectory above y=0y=0 . I realize that the title of this section is "codeless solution" and I'm about to show some code, but that's showbiz baby.
using Makie, CairoMakie

const y0 = 100.0
const v = 20.0
const g = 9.81

# keep θ in degrees
y(x; θ) = y0 + x*tan(π*θ/180) - g*sec(π*θ/180)^2*x^2/(2*v^2) 

xs = range(0.0, 120.0, step = 0.1)
ys = map(x -> y(x, θ = 20), xs)

fig = Figure()
ax = Axis(
    limits = (0, 120, 0, 110),
    title = "A firecracker trajectory", 
    xlabel = "horizontal distance from start (m)", 
    ylabel = "distance above ground (m)"
lines!(ax, xs, ys, color = :blue, label = L"\theta = 20^\circ")

Enter fullscreen mode Exit fullscreen mode

A parabola

We see a parabola, which is reassuring. But what's all this about "volume" of the region the projectile moves through? The question advises us that "It breaks into a large number of very small fragments, which move in every direction...", so let's try plotting many trajectories for various launch angles.

fig = Figure()
ax = Axis(
    limits = (0, 120, 0, 130),
    title = "Many firecracker trajectories", 
    xlabel = "horizontal distance from start (m)", 
    ylabel = "distance above ground (m)"

for θ in range(-90.0, 90.0, step = 3.0)
    xs = range(0.0, 120.0, step = 0.1)
    ys = map(x -> y(x, θ = θ), xs)
    lines!(ax, xs, ys)

Enter fullscreen mode Exit fullscreen mode

Many parabolas on one graph, they fill the space and appear to be "stuck" under some larger curve

Now the game is clear: the continuum of trajectories traces out a shape, so what we really need is somehow the "boundary" of all of these individual trajectories. If we knew that, we could determine the area traced out by all the trajectories, which is one step closer to the volume. But how to obtain such a thing?

To answer this, we have to think about what exactly is meant by "boundary" of the trajectories. Well, it's certainly some kind of function B(x)B(x) . But, given some xx , what is the value of B(x)B(x) ? If we examine the figure, the "boundary" of the trajectories at a given xx looks like the largest value attained by y(x)y(x) for any θ\theta . But that's exactly it:

B(x)=maxθ[π/2,π/2]y(x;θ) B(x) = \max_{\theta \in [-\pi/2, \pi/2]}{y(x; \theta) }

So for a given xx , we just need to compute the maximum of yy , treating it then as a function of θ\theta . This is accomplished by routine calculus,

0=yθ=xsecθ2gtanθsec2θv2x2    tanθ=v2gx 0 = \frac{\partial y}{\partial \theta} = x \sec \theta^2 - \frac{g \tan \theta\sec^2 \theta}{v^2} x^2 \implies \tan \theta = \frac{v^2}{g x}

Substituting this into y(x,θ)y(x, \theta) gives our boundary B(x)B(x) as
B(x)=y0+v22gg2v2x2 B(x) = y_0 + \frac{v^2}{2g} - \frac{g}{2 v^2} x^2

So (fascinatingly, in my opinion) the boundary of a set of parabolic trajectories is itself a parabola! The boundary is more usually called the envelope and is obtained for a parametric function by simply extremizing the parameter in question (in our case, θ\theta ) but I certainly didn't know that back in the day.

A neat observation: the quantity v22g\tfrac{v^2}{2g} is actually the height attained by a projectile shot straight up in the air. Therefore, setting ymax=v22gy_{\text{max}} = \tfrac{v^2}{2g} actually allows us to write

B(x)=y0+ymax(1x24ymax2) B(x) = y_0 + y_\text{max}\left(1 - \frac{x^2}{4y_\text{max}^2}\right)

Anyway, let's check out B(x)B(x) on our plot.

const ymax = v^2/(2*g)

B(x) = y0 + ymax*(1 - x^2/(4*ymax^2)) 

xs = range(0.0, 120.0, step = 0.1)
ys = B.(xs)
lines!(ax, xs, ys, 
    linewidth = 3.0, 
    color = :red, 
    label = L"\mathrm{Boundary: } B(x)")

Enter fullscreen mode Exit fullscreen mode

The same as the previous graph except now there is a thick red line on the boundary of all of the parabolas

Well there it is. So the area swept out by the trajectories is the area between the xx axis and the function B(x)B(x) . But how to get the volume? This is, to date, the absolute only single time I have used a volume of revolution calculation. You may have realized that the calculations we've done so far only hold for a "slice" of the volume. After all, the firecrackers are being ejected in all directions in space. But then, since each direction looks the same as any other, we can get the total volume from the area under B(x)B(x) by revolving B(x)B(x) around the yy axis using the method of cylinders. In brief, this simply states that the volume of a little chunk of cylinder is its circumference 2πr2 \pi r times its height hh times its width dr\text{d} r . Noting that B(2ymax1+y0/ymax)=0B(2y_\text{max} \sqrt{1 + y_0/y_\text{max}}) = 0 , we therefore integrate

Volume=02ymax1+y0/ymax2πxB(x)dx=2π(1+y0ymax)2ymax3, \text{Volume} = \int_{0}^{2y_\text{max} \sqrt{1 + y_0/y_\text{max}}} 2\pi x B(x) \,\text{d}x = 2 \pi \left(1 + \frac{y_0}{y_\text{max}} \right)^2 y_\text{max}^3 ,

a very slick result indeed. We can then obtain the final answer:

2*π*(1 + y0/ymax)^2 * ymax^3
Enter fullscreen mode Exit fullscreen mode
Enter fullscreen mode Exit fullscreen mode

The Julia Solution

To solve the problem in a "code-y" way, we'll essentially perform the above calculation, but numerically. First, we'll get B(x)B(x) by solving the optimization problem using Optim.jl. Then, we can integrate the result using QuadGK.jl.

using Optim
using QuadGK

boundary(x) = -1*optimize(θ -> -y(x, θ = θ), 0.0, 90.0).minimum
@time quadgk(x -> 2*π*x*boundary(x), 0, 2*ymax*sqrt(1 + y0/ymax))[1]
Enter fullscreen mode Exit fullscreen mode
  0.314286 seconds (286.15 k allocations: 16.542 MiB, 5.87% gc time, 99.95% compilation time)

Enter fullscreen mode Exit fullscreen mode

It's just that easy folks.


Even though this problem doesn't end up being that hard, it's been living in my head more or less rent-free for over ten years so it's nice to finally write it down. My Project Euler career has been over for a long time, but there's still some good stuff in there. Let me know if solving them gets you rich and famous.

Top comments (2)

3j5a profile image

how come you're neither famous nor rich :exploiding_head

great post!

gage profile image
Gage Bonner

Thank you!!