Iteration and convergence

In yesterday’s post, I had to find an angle, [\theta], for which

[\frac{1}{3} \theta = \sin \theta]

I wasn’t interested in the trivial solution, [\theta = 0], nor was I interested in the solution for which [\theta < 0]. I wanted the one that’s near [3 \pi/4], which turns out to be (to seven digits) 2.278863. How did I get it?

There are lots of ways to solve nonlinear equations like this. Newton’s method—or as I learned it, the Newton-Raphson method—is popular, but it does involve taking derivatives and then doing some algebra to set up the recurrence relation. And when there’s already a term that’s linear in the unknown variable, it’s really easy to set up a relation for fixed-point (sometimes called direct) iteration. Even though you usually need to work through more iterations than with Newton’s method, it’s often easier to do each iteration, especially when you’re tapping out the formula on a calculator.

In this problem, we can just move the 3 to the other side of the equation and iterate this way:

[\theta_{n+1} = 3 \sin \theta_n]

We start with an initial estimate, [\theta_0], take its sine, multiply by three, and that gives us the next estimate in the sequence.

Unfortunately, this doesn’t work. Even if you start with a initial estimate very close to the solution, say [\theta = 2.3], the fixed-point iteration sequence will never converge. Try it. You’ll see the values bouncing around the solution and getting further away with each cycle.

If we plot the two sides of the equation, their intersection will be the solution. We can also use this graph to track how the iteration proceeds.

Diverging iteration

As you can see, if we start with [\theta_0] that’s a little above the solution, we get a [\theta_1] that’s further below the solution than [\theta_0] was above it. Plugging that into the right-hand side gives us a [\theta_2] that’s back above the solution but even further away. And so on. So it looks like we’re SOL on fixed-point iteration for this problem.

Interestingly, if the problem were slightly different,

[\frac{1}{2} \theta = \sin \theta]

we would be able to use fixed-point iteration to get to the solution. Here, the recurrence relation is

[\theta_{n+1} = 2 \sin \theta_n]

and the sequence of iterations could be plotted like this,

Converging iteration

which is clearly heading in toward the solution.

The difference between the two problems boils down to the slope of the right-hand side term in the neighborhood of the solution. In our first problem, the slope at the solution is

[3 \cos 2.278863 = -1.951]

And in the second problem, the solution is [\theta = 1.895494] and the slope at the solution is

[2 \cos 1.895494 = -0.638]

The absolute value of the slope of the first is greater than one and the absolute value of the slope of the second is less than one, and that’s why the first diverges and the second converges. You can find the proof in the last two slides of these lecture notes for a class on numerical analysis.

So how did I find the solution to our original problem? I rearranged the problem as

[\sin^{-1} \frac{\theta}{3} = \theta]

from which I got the recurrence relation

[\theta_{n+1} = \pi - \sin^{-1} \frac{\theta_n}{3}]

If you’re wondering where the [\pi] term came from, that was necessary because the arcsine function on my calculator—like every calculator I’ve ever used—returns values in the first quadrant when the argument is positive. But our angle is in the second quadrant, so I had to take the supplement of what the arcsine returned.

I confess I didn’t bother taking the derivative of

[\sin^{-1} \frac{\theta_n}{3}]

to see if its absolute value was less than one when [\theta] was near [3 \pi/4]. I just assumed that if the slope of the sine was greater than one, the slope of its inverse would have to be less than one. So I entered a number near the solution and started iterating. It was easy to see after the first couple of cycles that the sequence was converging, so I kept at it until the answer stopped changing.

I don’t recommend this method of solving nonlinear equations in general. But I was watching the Standup Maths video that got me started on this rod rotation problem, and saw, starting around the 11:50 mark, that Matt was somewhat skeptical of Hugh’s solution for [\theta]. I didn’t want to go digging into SciPy to check the answer, so I opened PCalc and started iterating. My starting point was 2.2789, which was Hugh’s solution. It didn’t take long to get to 2.278863, which showed that Hugh’s answer was right to five digits.

I would have put this in yesterday’s post, but it was already long enough to stretch your patience. Now I have two posts that stretch your patience.

Pencil pushers

A couple of days ago, Matt Parker posted this video to his Standup Maths YouTube channel.

In it, he and Hugh Hunt, reader in engineering dynamics and vibration at Cambridge, stand a pencil on its end at the edge of a table, give the bottom end of it a horizontal whack, and see what happens. They also calculate what happens… sort of.

See, the thing is, the title of the video is “Will a falling pencil hit the table? We do the maths!” They certainly do some math, but not enough to answer the question in the title. Kind of disappointing. So let’s do a little more math and finish the problem for them. In the following, we’ll be referring to their work, so you really should watch the video before proceeding.

Here’s the setup. A pencil, which we are modeling as a uniform rod, is standing upright at the edge of a table. We’re going to apply a sharp, short duration horizontal force at the bottom of the rod and work out its subsequent motion. The sketch below shows the layout and the coordinate system I want to use in the solution.

Rod setup

The [g] represents the acceleration due to gravity. The rod is characterized by its length, [2a], and its mass, [m]. Because the mass is uniformly distributed along its length, the center of mass, [C], is halfway up the rod, and we can calculate its moment of inertia about [C] as

[I_C = \frac{1}{3} m a^2]

(If you’re looking at that formula and your first thought is “Why isn’t that a 12 in the denominator?” join the club. It took me a few minutes to remember that Hugh had defined [a] as half the length of the rod, and that’s why there’s a 3 instead of a 12. We’re used to seeing the formula written in terms of the full length

[I_C = \frac{1}{12} m L^2]

but since [L = 2a],

[I_C = \frac{1}{12} m (2a)^2 = \frac{1}{12} m (4a^2) = \frac{1}{3} m a^2]

so the 3 is correct.)

[A] is the end that’s initially at the top, and we’ll be tracking its position to find out if hits the table or flies past it.

[\hat{\:F}] is the impulse with which we strike the bottom of the rod. It’s shown as acting a little above the bottom, but I drew it that way so the arrow wouldn’t get lost in the top of the table. Impulse is defined as the integral of the force over the duration of its application. We’re going to assume the duration of the impulse is short, [\Delta t], and that it ends at time [t = 0]. Thus, the impulse is

[\hat{\:F} = \int_{-\Delta t}^0 F(t)\,dt]

The impulse-momentum equation is essentially just Newton’s second law integrated over time. For this situation, in the horizontal direction, we have

[\hat{\:F} = m v_0 - m v_{-}]

where [v_0] is the velocity of C at [t = 0] (just after the impulse), and [v_{-}] is the velocity of C at [t = -\Delta t] (just before the impulse). Because the rod starts out stationary, [v_{-} = 0] and

[\hat{\:F} = m v_0]

There’s a rotational analog to the linear impulse-momentum equation

[\hat{\:M}_C = \hat{\:F} a = I_C \,\omega_0 - I_C \,\omega_{-}]

where [\omega_0] is the angular velocity of the rod at [t = 0] and [\omega_{-}] is the angular velocity at [t = - \Delta t]. Because the rod is stationary before we hit it, [\omega_{-} = 0], so

[\hat{\:F} a = I_C \omega_0 = \frac{1}{3} m a^2 \omega_0]

There’s some simplification here that the video doesn’t mention. The full definition of angular impulse about [C] is

[\hat{\:M}_C = \int_{- \Delta t}^0 F(t) \, r_{\perp}(t) \, dt]

where [r_{\perp}(t)] is the perpendicular distance from the line of action of [\hat{\:F}] to [C]. It’s because the impulse is assumed to act over a very short length of time that we can take [r_{\perp}(t) = a] throughout the impulse. The idea is that impulse is so short, the rod doesn’t have a chance to rotate significantly. And because [a] is not a function of [t], we can bring it out of the integral, leaving us with [\hat{\:M}_C = \hat{\:F} a].

Combining our two impulse-momentum equations to eliminate [\hat{\:F}] gives us

[v_0 = \frac{1}{3} a \omega_0]

which is exactly what the video got.

Now let’s track the motions of [C] and [A]. After the impulse, the bottom of the rod will be just off the edge of the table, and the only force acting on the rod will be gravity (in the grand tradition of all elementary mechanics problems, we’re ignoring air resistance). So [C] will trace out a parabola, with

[x_C = v_0 \,t = \frac{1}{3} a \,\omega_0\, t]


[y_C = a - \frac{1}{2} g \,t^2]

being the parabola’s parametric equations.

With the only force acting through the rod’s center of mass, the angular velocity will remain constant at [\omega_0]. So

[\theta = \omega_0 \,t]

and we can use that to write the equations for the coordinates of [A]:

[x_A = x_C - a \sin \omega_0 t = \frac{1}{3} a\, \omega_0\, t - a \sin \omega_0 t] [y_A = y_C + a \cos \omega_0 t = a\,(1 + \cos \omega_0 t) - \frac{1}{2} g\, t^2]

Now, taking a hint from the video, we’re going to consider the conditions at which the tip of the rod, [A], is even, horizontally, with the edge of the table. In other words, when [x_A = 0].

[x_A = \frac{1}{3} a\, \omega_0\, t - a \sin \omega_0 t = 0]

That means

[\frac{1}{3} a\, \omega_0\, t = a \sin \omega_0 t]


[\frac{1}{3} \theta = \sin \theta]

This is the formula they got in the video, but they were looking only at the special case for which the tip of the pencil just barely hits the corner of the table. What we’ve shown is that the same formula works whenever the tip of the pencil passes a vertical line even with the edge of the table. And instead of stopping here, as the video does, we can use the solution to this formula to go further and answer the question in the video’s title: Under what conditions will the pencil hit the table?

The solution of the formula, to an unnecessarily precise seven digits, is [\theta = 2.278863], which is equivalent to about 130.57°. What we need to do is consider the value of [y_A] when [\theta] is equal to this value. There are three conditions of interest:

  1. If [y_A > 0] when [\theta = 2.278863], then the tip of the pencil will pass above the table and there will be no hit.
  2. If [y_A = 0] when [\theta = 2.278863], then the tip of the pencil will just graze the corner of the table. This is the boundary between hit and miss and is the condition they focused on in the video.
  3. If [y_A < 0] when [\theta = 2.278863], then the pencil hits the table.

Rod positions

Yes, [y_A] being less than zero for [\theta = 2.278863] is physically impossible because the pencil will hit the table before [\theta] can get to that value. But that’s OK. All we really need to know is the condition for which the pencil will hit the table. We can’t expect to track the motion of the pencil after it bounces off the table.

So let’s do some algebra on our expression for [y_A]. We can use the [\theta = \omega_0 t] relationship to eliminate [t] from the right-hand side:

[y_A = a\, (1 + \cos \theta) - \frac{1}{2}\frac{g}{\omega_0^2}\, \theta^2]

Dividing through by [a], we get a formula with nondimensional terms, which will allow us to do some numerical work.

[\frac{y_A}{a} = (1 + \cos \theta) - \frac{1}{2}\frac{g}{a\, \omega_0^2}\, \theta^2]

Substituting in [\theta = 2.278863] and rearranging, we get these three conditions:

  1. [y_A > 0] when [\frac{g}{a\, \omega_0^2} < 0.134650] (misses table)
  2. [y_A = 0] when [\frac{g}{a\, \omega_0^2} = 0.134650] (grazes table)
  3. [y_A < 0] when [\frac{g}{a\, \omega_0^2} > 0.134650] (hits table)

So for a given rod length and a given gravitational field, it’s the angular velocity imparted to the rod by the impulse that determines whether the rod strikes the table or not. If the angular velocity is above this limit,

[\omega_0 = \sqrt{\frac{g}{0.134650\, a}}]

the rod will not hit the table. If the angular velocity is below this value, it will.

We can now work out what this means in terms of the magnitude of the impulse. Recall that

[\hat{\:F} = \frac{1}{3} m\, a\, \omega_0]

so we can solve for the value of [\hat{\:F}] that is the boundary between hit and miss:

[\hat{\:F} = 0.908396\, m\, \sqrt{g\, a}]

Impulses higher than this will be misses, impulses lower will be hits.

If you stop for a minute and think about this result, you should be able to convince yourself that it makes sense. Clearly, a larger mass will require a stronger impulse to get the rod to fly off the table. And a higher gravitational field would mean the rod drops more quickly, also requiring a stronger impulse to get it to miss the table. The relationship between [\hat{\:F}] and [a] is perhaps not as obvious, but think of it this way: The longer the rod, the more slowly it rotates for a given impulse; therefore, the only way to get it to swing around to the required angle of [\theta = 2.278863] before striking the table is to hit it harder.

Now let’s run some numbers for what we see in the video. We’ll assume the pencil is of length 18 cm. With [g = 9.81\:\mathrm{m/s^2}] and [a = 0.09\:\mathrm{m}], we get [\omega_0 = 28.45\:\mathrm{rad/s}] (about 270 rpm) as the boundary between hit and miss. That corresponds to

[v_0 = \frac{1}{3}a\, \omega_0 = 0.85\:\mathrm{m/s}]

as the horizontal component of velocity at [C]. The velocity of the base of the pencil just after the impulse needs to be

[v_0 + a\, \omega_0 = \frac{4}{3} a\, \omega_0 = 3.4\:\mathrm{m/s}]

These values seem pretty easy to achieve, which is probably why most of the hits in the video sent their pencils flying off without hitting the table.

Larry Tesler

You probably read a few articles this past week commemorating Larry Tesler, the father of cut/copy-paste and the most prominent advocate of modeless computing. It’s good to read what other have to say about him, but if you have the time you should read some of his own writing.1

His ACM article, “A Personal History of Modeless Text Editing and Cut/Copy-Paste,” is probably the best distillation of the ideas he’s known for and the connection between them. The main thread is a fun historical narrative, but don’t skip over the sidebars in the article. In “How Modes Degrade Usability,” he ties modes to the once-predominant verb/object method of telling a computer what to do. In verb/object computing, you give the computer an action command and follow that with the item(s) on which to perform the action. On a computer of today, that would be like selecting Copy from the Edit menu and then selecting the text to copy.

That may sound insane, but one of the things that made the Mac initially puzzling to people who were heavy users of PCs and mainframes back in the mid-80s was that it flipped verb/object around. I remember having to explain to several people that the main thing to remember on a Mac was that you never go to the menu bar for a command until you’ve selected the thing that the command would act on. The verb/object ordering was hard to shake for them because it seemed very English-like and conversational. Those of us who use the command line (or Vim2) still do a lot of verb/object computing, but nowadays we do it within a modeless windowing environment. Modal computing tends to be the exception rather than the rule.

Tesler’s ideas about how computing should be done has become so pervasive I imagine it’s almost inconceivable to younger readers that there was ever a different way. When Steve Jobs stepped down as CEO of Apple, I wrote “In a sense, you’re using a Steve Jobs product whether it has an Apple logo or not.” Much the same can be said of Larry Tesler.

  1. I really like his annotated version of the PUB manual, but I have a thing for markup languages. 

  2. Kakoune is a modal text editor that turns verb/object around. By highlighting selections (like Vim’s visual mode) all the time, it allows you to make selections first, see what they are, and then issue the commands that act on them. It’s an interesting exception to the connection Tesler draws between modal and verb/object ordering. 

A name display bug in Files

I’ve written in the past about how I’ve had to use tags to work around a misfeature of Files: that it doesn’t display full file names, including the extensions. Today I learned that there’s a bug in that misfeature.

I was working from my iPad on some files that were stored in iCloud Drive. Here’s Files showing the contents of the folder I was working in:

Files inconsistent in showing extensions

The file with the blue dot has a .tex extension. The extension isn’t shown, which is the normal (and frustrating) way Files behaves. But the file with the green dot has a .py extension and its extension is shown. What’s going on?

My first concern was that the green file somehow had two extensions—that its full name was—so I opened a Terminal in Textastic,1 navigated to the appropriate iCloud folder, and ran ls to see the real, full names of both files.

Terminal session showing full file names

As you can see, there was no double extension.2 The problem wasn’t that I had somehow misnamed the Python file, it was that Files has a bug in the name display. It happens to be a bug I approve of, and one that I wish were consistent across all file types at all times, but it’s a bug nonetheless.

It says something about the state of software quality at Apple that it can’t even get its design errors right.

  1. Normally, I’d do this in Prompt, but I was in the middle of editing one of the files in Textastic when I noticed the extension anomaly in Files, so it was quicker to start a terminal session there. I really like Textastic. 

  2. If you’re looking at the path in the terminal session and wondering why an iCloud folder isn’t in ~/Library/Mobile Documents/com~apple~CloudDocs, its because I made a symbolic link to the projects folder in iCloud from a projects directory in my home directory. Saves time when using the cd command, whether I’m working directly on my Mac or remotely through terminal emulation.