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.