Saturday, August 3, 2013

The 10° Solution

Last week we looked at how to get analytic expressions for the sine and cosine of angles between 0 and 90 °, provided that the angle was an integer multiple of 3°.

For aesthetic reasons, if nothing else, it would be nice to know have a table of exact expressions for sines and cosines for every degree, not just every 3°. In principle we could do this, following the tricks we did last time, by noting that

      cos 90° = 0  ,  (1)

and then expand

      cos 90 x = Σn=090 an cosn x = 0  ,  (2)

which gives us a 90th degree polynomial in cos x, one of whose solutions is cos 1°.

Unfortunately (2) is just a tad difficult to solve. But there are other ways to do this. If, for example, we knew cos 10°, then we could get cos 1° and sin 1° from the addition formulas:

     cos 1° = cos 10° cos 9° + sin 10° sin 9° ,  (3)

     sin 1° = sin 10° cos 9° - sin 9° cos 10°  .  (4)

So how to we find cos 10°, or, in radians, cos π/18 ? Since we know cos π/2 = 0, it follows that π/18 will be one of the solutions of the equation

      cos 9 x = 0  .  (5)

All we have to do is find that solution.

We start by expanding (5), writing

      cos 9 x = cos x (3 - 4 cos2 x) ( 3 - 36 cos2 x + 96 cos4 x - 64 cos6 x )  .  (6) *

Let's consider the possible solutions we can get for (6). We know that cos x is zero at x = π/2, and at multiples of π away from this point, so our solutions are going to be of the form

      9 x = ½ π + N π  ,  (7)

which gives 9 unique solutions:

      x = π/18, π/6, 5 π/18, 7 π/18, π/2, 17 π/18, 5 π/6, 13 π/18, and 11 π/18 .  (8)

Since cos(π-x) = - cos x, the cosines at the last four points here will be just the negatives of the cosines of the first four points, which is not surprising, since (6) obviously has a lot of ± roots.

Given that the cosine decreases as x goes from 0 to ½π, cos 10° will be the largest root of (6).

The first term of (6) vanishes when x = ½π, while the second term vanishes when x = π/6 or 5π/6. So all we have to worry about is the third term. If we set

      y = cos2 x ,  (9)

then we only need find the roots of the third degree polynomial:

      P(y) = 64 y3 - 96 y2 + 36 y - 3  .  (10)

There is a standard method for doing this, which dates back almost 500 years. I'm going to work it out, once, because, well, this is Stupid Math Tricks, and I've never written it down before.

The first step is to get rid of the quadratic term in (10). To do this, set y = z + δ and collect powers of z:

      P(z + δ) = 64 z3 + 96 (2 δ - 1) z2 + 12 (16 δ2 - 16 δ + 4) z + (64 δ3 - 96 δ2 + 36 δ - 3)  .  (11)

If we set δ = ½ we get

      P(z + ½) = 64 z3 - 12 z - 1  .  (12)

(You can check this is right by noting that P(0) = -3, P(1) = 1, P''(0) = 36, and P'''(0) = 384 in both (11) and (12).)

Now write

      z = u + v  .  (13)

We'll collect the terms in a somewhat non-obvious way: write

      z3 = u3 + v3 + 3 u v (u + v)  ,  (14)

and collect terms in (u+v):

      P(u + v + ½) = 64 u3 + 64 v3 + 12 (16 u v - 1) (u + v) - 1  .  (15)

Get rid of the third term in (15) by setting

      v = 1/(16 u)  .  (16)

Then

      P(u + 1/(16 u) + ½) = 64 u3 + 1/(64 u3) - 1  .  (17)

Solving for P(y) = 0 is therefore equivalent to solving the equation

      4096 u6 - 64 u3 + 1 = 0  .  (18)

This is a quadratic equation in u3. There are no real solutions, but it does have the complex solutions:

      u3 = (1 ± i √3)/128 = (1/64) ( cos π/3 ± i sin π/3 ) = (1/64) exp ±i π/3  .  (19)

This gives us

      u = (1/4) exp ±i π/9  .  (20)

Without loss of generality we can take the + sign in (20), in which case we find

      u = (1/4) exp i π/9

      v = (1/4) exp -i π/9  .  (21)

Then

      u + v = (1/4) (exp i π/9 + exp -i π/9) = ½ cos π/9  .  (22)

Or, to get back to the solution of (10),

      y = z + ½ = u + v + ½ = ½( 1 + cos π/9)  .  (23)

Using the half-angle formula for cosines, we find

      y = cos2 π/18  .  (24)

And so one of the solutions of (5) is TaDa!

      cos x = √y = cos π/18  .  (25)

But, of course we know that one of the solutions of (6) must be cos π/18, so we've come a long way to prove:

      cos π/18 = cos π/18  .  (26)

THAT's a big disappointment. While cos π/18 is an algebraic irrational, i.e. not transcendental, it cannot be expressed in any simple formula such as those we found last time.

You can, however, with enough patience get a solution. Using Newton's method, we can find a root of any differentiable function f(x) by setting x0 sufficiently close to the root and repeatedly calculating

      xn+1 = xn - f(xn)/f'(xn)  .  (27)

Let's start with the simplest polynomial we have, (12), which we'll write as

      F(z) = 64 z3 - 12 z - 1  .  (28)

The corresponding Newton function is

      G(z) = z - F(z)/F'(z) = (128 z3 + 1)/[12 (16 z2 + 1) ]  .  (28)

Now we know that z = y - ½ = cos2 π/18 - ½, and that cos π/18 is near one, so a good starting value for z is ½. That let's us right this little awk script. I picked awk because it's even available in Windows. I'm using bash to write the script, so you'll have to modify that for Windows systems. It will work on Macs and Linux boxes:

#! /bin/bash

Z=0.5

# A proper script would turn off when convergence is achieved,
#  and not before.  Obviously this is not a proper script

# This is set near the maximum precision of my computer. YMMV

for i in `seq 1 5`
do
    Z=`echo $Z | awk '{printf "%21.15f", (128*$1^3+1)/(12*(16*$1^2-1))}'`
    echo $Z
done

Y=`echo $Z | awk '{printf "%21.15f", $1+0.5}'`

echo "Y = " $Y

COS=`echo $Y | awk '{printf "%21.15f", sqrt($1)}'`

echo "cos pi/18 = " $COS

Which has the output:

0.472222222222222
0.469862891737892
0.469846311209168
0.469846310392954
0.469846310392954
Y =  0.969846310392954
cos pi/18 =  0.984807753012208

Now to go back to last week's barbarian captor scenario, since you'll be doing this with a stick on dirt it's going to take a while, but you can tell him that, in principle, you can get the exact answer to arbitrary precision.


* You can get this formula by

  • Carefully working out the expansion of cos 9x using the addition formulas, then factoring the resulting polynomial.
  • Or realize that cos (π - x) = - cos x means that any expansion of cos 9x must be an odd polynomial in cos x, and that x = π/6 is one of the solutions of (5), so that ½ √3 must be one of the solutions, and derive the remaining coefficients by careful consideration of how cos 9x behaves, e.g., both sides of (5) must be equal to 1 when x = 0. Or,
  • Or do what I did: write out (5) with unknowns for the variables, use a spreadsheet, or Perl or a Fortran program to write down cos 9x for one hundred or so values of x, and use gnuplot's fitting function to determine the values of the coefficients. [Go back to equation (6).]

† Actually, since exp 2 N i π = 1 for any integer N, there are also two other unique solutions, u = exp ±i 7π/9 and u = exp ±i 5π/9. We'll leave it as an exercise for the reader to follow through with these solutions. (Hint: the final solution will involve the cosines o 10°, 50°, and 70°.) [Go back to equation (20).]