Math for Programmers 3D graphics, machine learning, and simulations with Python

Testing a projectile simulation 429 From this plot, we can guess what the optimal value is. At a launch angle of 45°, the landing position is maximized at a little over 40 meters from the launch position. In this case, 45° turns out to be the exact value of the angle that maximizes the landing position. In the next section, we’ll use calculus to confirm this maximum value with- out having to do any simulation. 12.1.4 Exercises Exercise 12.1 How far does the cannonball go when fired at an angle of 50° from an initial height of zero? How about if it is fired at an angle of 130°? Solution At 50°, the cannonball goes about 40.1 meters in the positive direc- tion, while at 130°, it goes 40.1 meters in the negative direction: >>> landing_position(trajectory(50)) 40.10994684444007 >>> landing_position(trajectory(130)) -40.10994684444007 This is because 130° from the positive x-axis is the same as 50° from the negative x-axis. Exercise 12.2—Mini Project Enhance the plot_trajectories function to draw a large dot on the trajectory graph at each passing second so we can see the passing of time on the plot. Solution Here are the updates to the function. It looks for the index of the nearest time after each whole second and makes a scatter plot of (x, z) values at each of these indices: def plot_trajectories(*trajs,show_seconds=False): for traj in trajs: xs, zs = traj[1], traj[2] plt.plot(xs,zs) if show_seconds: second_indices = [] second = 0 for i,t in enumerate(traj[0]): if t>= second: second_indices.append(i) second += 1 plt.scatter([xs[i] for i in second_indices], [zs[i] for i in second_indices]) ...

430 CHAPTER 12 Optimizing a physical system (continued) As a result, you can picture the elapsed time for each of the trajectories you plot; for example, plot_trajectories( trajectory(20), trajectory(45), trajectory(60), trajectory(80), show_seconds=True) 20.0 80° 60° 17.5 0 Height 15.0 12.5 10.0 20° 45° 40 7.5 20 5.0 Distance 2.5 0.0 10 30 Plots of four trajectories with dots showing their positions at each whole number of seconds. Exercise 12.3 Make a scatter plot of hang time versus angle for angles between 0 and 180°. Which launch angle produces the maximum hang time? Solution test_angles = range(0,181,5) hang_times = [hang_time(trajectory(theta)) for theta in test_angles] plt.scatter(test_angles, hang_times) Hang time (seconds) 4.0 A plot of the hang time of 3.5 the cannonball as a function 3.0 of the launch angle 2.5 2.0 25 50 75 100 125 150 175 1.5 θ 1.0 0.5 0.0 0

Testing a projectile simulation 431 It appears that a launch angle of roughly 90° yields the longest hang time of just about 4 seconds. This makes sense because  = 90° yields the initial velocity with the largest vertical component. Exercise 12.4—Mini Project Write a function plot_trajectory_metric that plots the result of any metric we want over a given set of theta ( ) values. For instance, plot_trajectory_metric(landing_position,[10,20,30]) makes a scatter plot of landing positions versus launch angle for the launch angles 10°, 20°, and 30°. As a bonus, pass the keyword arguments from plot_trajectory_metric to the internal calls of the trajectory function, so you can rerun the test with a different simulation parameter. For instance, this code makes the same plot but simulated with a 10-meter initial launch height: plot_trajectory_metric(landing_position,[10,20,30], height=10) Solution def plot_trajectory_metric(metric,thetas,**settings): plt.scatter(thetas, [metric(trajectory(theta,**settings)) for theta in thetas]) We can make the plot from the previous exercise by running the following: plot_trajectory_metric(hang_time, range(0,181,5)) Exercise 12.5—Mini Project What is the approximate launch angle that yields the greatest range for the cannonball with a 10-meter initial launch height? Solution Using the plot_trajectory_metric function from the preceding mini-project, we can simply run plot_trajectory_metric(landing_position,range(0,90,5), height=10)

432 CHAPTER 12 Optimizing a physical system (continued) 50 40 Range 30 20 10 0 20 40 60 80 Launch angle A plot of range of the cannonball vs. launch angle with a 10 meter launch height The optimal launch angle from a height of 10 meters is about 40°. 12.2 Calculating the optimal range Using calculus, we can compute the maximum range for the cannon, as well as the launch angle that produces it. This actually takes two separate applications of calculus. First, we need to come up with an exact function that tells us the range r as a function of the launch angle  . As a warning, this will take quite a bit of algebra. I’ll carefully walk you through all the steps, so don’t worry if you get lost; you’ll be able to jump ahead to the final form of the function r( ) and continue reading. Then I show you a trick using derivatives to find the maximum value of this func- tion r( ), and the angle  that produces it. Namely, a value of  that makes the deriv- ative r'( ) equal zero is also the value of  that yields the maximum value of r( ). It might not be immediately obvious why this works, but it will become clear once we examine the graph of r( ) and study its changing slope. 12.2.1 Finding the projectile range as a function of the launch angle The horizontal distance traveled by the cannonball is actually pretty simple to calcu- late. The x component of the velocity vx is constant for its entire flight. For a flight of total time Δt, the projectile travels a total distance of r = vx · Δt. The challenge is find- ing the exact value of that elapsed time Δt. That time, in turn, depends on the z position of the projectile over time, which is a function z(t). Assuming the cannonball is launched from an initial height of zero, the first time that z(t) = 0 is when it’s launched at t = 0. The second time is the elapsed time we’re looking for. Figure 12.9 shows the graph of z(t) from the simulation with 

Calculating the optimal range 433 = 45°. Note that its shape looks a lot like the trajectory, but the horizontal axis (t) now represents time. trj = trajectory(45) ts, zs = trj[0], trj[2] plt.plot(ts,zs) 10 z(t) 8 6 z 4 2 Landing time Δt ≈ 2.9s z(t) = 0 Launching time 0 z(0) = 0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 t Figure 12.9 A plot of z(t) for the projectile showing the launching and landing times where z = 0. We can see from the graph that the elapsed time is about 2.9 seconds. We know z''(t) = g = –9.81, which is the acceleration due to gravity. We also know the initial z velocity, z'(0) = |v| · sin( ), and the initial z position, z(0) = 0. To recover the position function z(t), we need to integrate the acceleration z''(t) twice. The first inte- gral gives us velocity: t z (t) = z (0) + gdτ = |v| · sin(θ) + gt 0 The second integral gives us position: z(t) = z(0) + t t |v| · sin(θ) + gτ dτ = |v| · sin(θ) · t + g t2 0 2 z (τ ) dτ = 0 We can confirm that this formula matches the simulation by plotting it (figure 12.10). It is nearly indistinguishable from the simulation. def z(t): A direct translation of the result of return 20*sin(45*pi/180)*t + (-9.81/2)*t**2 the integral, z(t), into Python code plot_function(z,0,2.9)

434 CHAPTER 12 Optimizing a physical systemz 10 8 6 4 2 Figure 12.10 Plotting the 0 exact function z(t) on top of the simulated values 0.0 0.5 1.0 1.5 2.0 2.5 3.0 t For notational simplicity, let’s write the initial velocity |v| · sin( ) as vz so that z(t) = vzt + gt2/2. We want to find the value of t that makes z(t) = 0, which is the total hang time for the cannonball. You may remember how to find that value from high school alge- bra, but if not, let me remind you quickly. If you want to know what value of t solves an equation at 2 + bt + c = 0, all you have to do is plug the values a, b, and c into the qua- dratic formula: √ t = −b ± b2 − 4ac 2a An equation like at2 + bt + c = 0 can be satisfied twice; both times when our projectile hits z = 0. The symbol ± is shorthand to let you know that using a + or – at this point in the equation gives you two different (but valid) answers. In the case of solving z(t) = vzt + gt 2/2 = 0, we have a = g/2, b = vz and c = 0. Plugging into the formula, we find t = −vz ± vx2 = −vz ± vz g g Treating the ± symbol as a + (plus), the result is t = (–vz + vz)/g = 0. This says that z = 0 when t = 0, which is a good sanity check; it confirms that the cannonball starts at z = 0. The interesting solution is when we treat ± as a – (minus). In this case, the result is t = (–vz – vz)/g = –2vz/g. Let’s confirm the result makes sense. With an initial speed of 20 meters per second and a launch angle of 45° as we used in the simulation, the initial z velocity, vz, is –2 · (20 · sin(45°))/–9.81 ≈ 2.88. This closely matches the result of 2.9 seconds that we read from the graph.

Calculating the optimal range 435 This gives us confidence in calculating the hang time Δt as Δt = –2vz/g or Δt = –2|v| sin( )/g. Because the range is r = vx · Δt = |v|cos( ) · Δt, the full expression for the range r as a function of the launch angle  is r(θ) = − 2|v|2 sin(θ) cos(θ) g We can plot this side by side with the simulated landing positions at various angles as in figure 12.11 and see that it agrees. def r(theta): return (-2*20*20/-9.81)*sin(theta*pi/180)*cos(theta*pi/180) plot_function(r,0,90) Landing position 40 r(θ) Figure 12.11 Our calculation of 35 projectile range as a function of 30 Simulations the launch angle r( ), which 25 matches our simulated landing 20 20 40 60 positions 15 θ 10 80 5 0 0 Having a function r( ) is a big advantage over repeatedly running the simulator. First of all, it tells us the range of the cannon at every launch angle, not just a handful of angles that we simulated. Second, it is much less computationally expensive to evalu- ate this one function than to run hundreds of iterations of Euler’s method. For more complicated simulations, this could make a big difference. Additionally, this function gives us the exact result rather than an approximation. The final benefit, which we’ll make use of next, is that the function r( ) is smooth, so we can take its derivatives. This gives us an understanding of how the range of the projectile changes with respect to the launch angle. 12.2.2 Solving for the maximum range Looking at the graph of r( ) in figure 12.12, we can set our expectations for what the derivative r'( ) will look like. As we increase the launch angle from zero, the range increases as well for a while but at a decreasing rate. Eventually, increasing the launch angle begins to decrease the range.

436 CHAPTER 12 Optimizing a physical system The key observation to make is that while r'( ) is positive, the range is increasing with respect to  . Then the derivative r'( ) crosses below zero, and the range decreases from there. It is precisely at this angle (where the derivative is zero) that the function r( ) achieves its maximum value. You can visualize this by seeing that the graph of r( ) in figure 12.12 hits its maximum when the slope of the graph is zero. 40r rʹ(θ) = 0 r (θ) 35 rʹ(θ) < 0 20 40 60 80 Range decreases 30 θ as angle increases 25 20 rʹ(θ) > 0 Range increases 15 as angle increases 10 5 0 0 Figure 12.12 The graph of r( ) hits its maximum when the derivative is zero and, therefore, the slope of the graph is zero. We should be able to take the derivative of r( ) symbolically, find where it equals zero between 0° and 90°, and this should agree with the rough maximum value of 45°. Remember that the formula for r is r(θ) = − 2|v|2 sin(θ) cos(θ) g Because –2|v|2/g is constant with respect to  , the only hard work is using the product rule on sin( )cos( ). The result is r (θ) = 2|v|2 cos2(θ) − sin2(θ) g Notice that I factored out the minus sign. If you haven’t seen this notation before, sin2( ) means (sin( ))2. The value of the derivative r'( ) is zero when the expression sin2( ) – cos2( ) is zero (in other words, we can ignore the constants). There are a few ways to figure out where this expression is zero, but a particularly nice one is to use the trigonometric identity, cos(2 ) = cos2( ) – sin2( ), which reduces our prob- lem even further. Now we need to figure out where cos(2 ) = 0.