Math Notes 26

On Newton’s Method

Newton’s method allows finding the ‘root’ of a function (where it is zero), given conditions expressed in the wikipedia article no doubt.

More than this however, it can be used to find the inverse of a function. The common example is $\sqrt{2}$. If you know $x$ where $x^2-2=0$, then you know $\sqrt{2}$.

The rough gist of the algorithm is as follows.

1. Take a function like $f(x) = x^2-2$
2. Make an arbitrary guess on the value of $x$ where $f(x)=0$, $g_0$
3. Evaluate $f(x)$ at $g_0$ and find the line tangent to $f(x)$ at $g_0$
4. Find the x-intercept of the tangent line. This x-intercept value is $g_1$
5. Repeat steps 2…5 (incrementing $n$ throughout) until a satisfactory precision is obtained.

Algebraically, we start with the function and its derivative:

And make a guess $g_0$ and create the tangent line at that guess location on $f(x)$

find $b$ by evaluating at $g_0$:

Which finishes the line equation as:

Next, find the x-coordinate by evaluating at $y=0$

As stated above, this $x$ value is actually $g_1$, or generally

Going back the case of $\sqrt{2}$, we substitute in the function

With $g_n=5$ we see the first four approximations as:

Generally…

This process is easily generalized to finding inverses instead of roots — just as in the $\sqrt{2}$ example.

$f(x)$ was previously defined as $x^2-2$, but more generically is described as $f(x)=x^2$. Then, inside the recursive equation, the $-2$ is placed outside the function

The above $g_{n+1}$ approximates $f^{-1}(2)$. the $2$ is replaced with some constant $c$ to find $f^{-1}(c)$:

The generic case for $\sqrt{c}$ then is:

Or

One way to loosely verify this is to: consider that when $n$ approaches $\infty$, the differences between $g_n$ and $g_{n+1}$ will grow increasingly smaller (assuming the limit does indeed converge) and indeed will be come so close that we can consider them to be equal. So, we set $g_{n+1} = g_n = x$ and solve:

Other inverses

This method can be extended to other functions. Some examples:

Approximate Functions

By expanding the recursion a few steps, we can create an obnoxious approximation for $\ln(x)$:

And one for $\sqrt{x}$:

In a way, this approximation for $\sqrt{x}$ is quite reasonably constructed. It has plenty of detail around the origin, but as it goes further the higher-order polynomial components divide out and all that is left is the slope asymptote.

These examples were generated with the Python library SymPy. (I got exhausted of writing it by hand around third recursion)

from sympy import symbols, simplify, print_latex

x = symbols("x")

expr = 1
for i in range(0, 3):
    expr = 0.5 * expr + x / (2 * expr)

print_latex(simplify(expr))

Newton’s Method: π

Suppose $f(x)=\sin(x)$, then $f(\pi) = 0$. Using Newton’s method, we arrive at the recursion:

Which gives π accurate to 100 digits in only 5 iterations.

from sympy import N, tan, pi

expr = 3
for i in range(0, 5):
    expr = expr - tan(expr)

print(N(expr, 100))
print(N(pi - expr, 100))

The behavior of $g_{n+1} = g_n - \tan(g_n)$ depends highly the initial value $g_0$. This makes some intuitive sense because $\sin(x)$ has many different roots so the sequence would likely converge to the nearest root to the starting point.

This behavior is clear by plotting $f(x) = x - \tan(x)$:

Around each multiple of pi, the graph flattens out. to the right of each flat point, the value of the function is decreases slightly. To the left, the value increases slightly. When the function is interpreted as the mapping from the previous value $x$ to the next value $f(x)$, the slight decreases and increases near each multiple of π makes the recursion stably approach the nearest multiple of π. Further away from the multiples of π, the slope becomes steep enough to ‘push’ the next step of the recursion away from the nearest multiple of pi. The steep region is not well-behaved. for these notes, ‘well-behaved’ means that the recursion $g_{n+1} = g_n - \tan(g_n)$ has a value which approaches the multiple of π nearest to $g_0$.

Specifically, let’s look at $x=0\pi$. The well-behaved region will extend an equal distance positive and negative from zero (because the function is symmetric to an extent around the multiples of π) Therefore, we can solve for this distance by finding $g_n = g_{n+1}$, or $x=-(x-tan(x))$.

Or

Where $g(x) = 0$

Plotting this function, we see the zero points are around $\approx\pm1.166$.

I can’t think of any algebra to solve $g(x) = 0$. Instead, we will use Newton’s method to find some ‘radius of stability’ $S$ such that $g(S)=0$

(there’s some humor in here somewhere; I just know it!)

Starting with $g_0 = 1.16$, we see it level out to $S\approx1.16556118$. I do not believe this constant has any simpler symbolic forms.

from sympy import N, tan, sec

expr = 1.16
for i in range(0, 15):
    expr = expr - (tan(expr) - 2 * expr) / (sec(expr) * tan(expr) - 2)
    # evaluate algebra to numeric value each step to prevent
    # algebraic explosion. Use 200 digits to avoid precision loss
    expr = N(expr, 200)
print(N(expr, 15))

So, around $0\pi$ this well-behaved region is described as

And in general is:

Where:

When $g_n$ is outside the well-behaved region, $g_{n+1}$ will skip to be nearby an entirely different multiple of pi. That makes me curious: is there some value of $g_{0}$ such that all subsequent $g_n$ values will keep skipping away and never converge?

Well, there’s the solution $g_{0} = S$ where the iterations bounce back and forth, but that’s more of a bounce than a skip. In theory, there exists some value $g_n$ such that $g_{n+1}=g_{n}+\pi$. A $g_n$ of this description will never converge.

Of course, this same logic applies for the situation where $g_{n+1} = g_{n} \pm q\pi$ where $q\in\mathbb{Z}$ and $q\ne 0$. So, generally:

will diverge. For example:

I was quite interested in describing the nature of this non-well-behaved region. However, the more I looked, the more detail I saw. Slight changes in input create vast differences in output, nudging the inputs even slightly creates entirely different outputs. This detail seems to persist no matter how precise the inputs are.

This all stinks of a fractal.

After discovering that the details were infinite, I stopped trying to directly define the results of the function. However, one last detail I noticed which I have not yet explained is the converging values around $-\pi/2$ (and all other midpoints between adjacent well-behaved areas)

As seen in the graph, the convergent values around $-\pi/2$ generally tend to be the multiple of π less than a certain hyperbola. Fascinating.

I recommend viewing this graph in Desmos yourself. I’ve included the code to recreate it below (make sure to make the range on the y-axis quite large).

f\left(x\right)=x-\tan\left(x\right)

x=-\frac{\pi}{2}

y=f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(x\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)

\frac{1}{\left(x+\frac{\pi}{2}\right)}

\operatorname{floor}\left(\frac{\frac{1}{\left(x+\frac{\pi}{2}\right)}}{\pi}\right)\pi

Perhaps this is quite similar to 3Blue1Brown: Newton’s Fractal. It has been a while since I have seen that video.

Proof that π is rational (joke)

from sympy import N, pi, floor
best = (10000, 1)

for i in range(1, 1_000_000):
    val = N(pi * i, 100)
    lz = val - floor(val)
    print(i, lz, val)
    lz = lz.evalf()
    if lz < best[0]:
        best = (lz, i)

print(best)

If I were to trust google’s calculator… then I’d say π is rational.

On A Spigot Algorithm for the Digits of π (Rabinowitz and Wagon)

Two equivalent expressions of a standard decimal number:

A mixed-radix number is somewhat like having a different base for each digit. The primary objective of this paper is stating mixed-radix representation’s of $e$ and $\pi$, and converting them to more useful bases.

Example of mixed-radix number:

Which expands to

When $a_i=1$ for all $i$, this is the well-known power series for $e$.

Euler’s number spigot

A. H. J. Sale doi 10.1093/comjnl/11.2.229.

Well-known $e$ power series:

Given in mixed-radix form:

Also can be stated as a recursion:

The paper then multiplies the fractional part of the expression by $\frac{10}{10}$, and distributes the numerator into the recursion:

The goal of this manipulation is to extract base-10 digits from this mixed-radix form, one digit at a time. We first take this recursion to some specific depth. $n=6$ in this example.

From right to left, we consecutively reduce the inner values mod the outer divisor, moving out the quotient in the process. An example is more clear:

This method of distributing the divisor, diving with remainder, and factoring is repeated:

This shows clearly that the 2nd digit is $7$. To get further digits, we again multiply the mixed-radix part by $10/10$.

And repeat the process:

And again:

And again:

Which now produces an incorrect result.

Why? Well, it is reasonable to assume that, because we limited the recursion to a finite number of levels, the remaining levels begin to become significant as the powers of $10$ rise.

Looking at the last reasonable result, we see at the end:

Let’s call $(\ldots)$ to be $r$ , we see

Which expands to

Multiplying in the $10^3$ creates a value less than $1$, and so it does not effect the calculations.

However, when the power of ten is increased to $10^4$, we see:

Which is greater than one and will likely cause a deleterious effect when it is missing from the calculations. So, the question is: how many terms does the recursion need to be expanded to to obtain $d$ digits? Instead of thinking about the algorithm as algebraic manipulations on a recursive equation, I find it easier to think of as an array of numbers. The first step of the above example becomes:

Where:

Inside the array, there are 5 slots (columns). Each slot represents a level of recursion. More slots => more digits extractable.

Let $d$ represent the number of digits being extracted, and $s$ represent the number of slots. A calculation using $s$ slots will yield $d$ digits when the following is true:

Infinite sums are cumbersome so instead a simpler form is desired. So long as the value of the sum overestimated, a simpler form will still work.

The following form is valid according to Desmos and WolframAlpha, so I’ll trust those sites on that.

Substituting this simpler expression in:

Solve for digits:

Plotting this, we see

The linear approximation of

Fits well under the entire domain. This approximation becomes poorer as more slots are added. When $s=70$, $d \le \approx100$

A polynomial or piecewise approximation could provide improvement. The original paper provides a different approximation entirely. I have not reviewed that approximation.

When in one of the previous examples the result of $2.7180\ldots$ was obtained, only 5 slots were used, and so only two digits after the decimal place should have been expected with confidence to be correct.

Back again to thinking of the problem algorithmically:

Where:

To find the second digit, update the table such that:

• $v_i=10r_i$ for all $i$
• $q_i$ is cleared for all $i$

Then apply the same rules, stopping after $i=1$:

And so digit #2 is $1$: $2.71\ldots$.

An implementation in Rust is provided.

fn main() {
    for i in 1..20 {
        println!("{:?}", e(i));
    }
}

fn e(d: u32) -> Vec<u8> {
    let mut digits = vec![0; d as usize];
    let slots = d as usize + 3;
    let mut tmp = vec![10; slots];
    for stop in 0..d as usize {
        let mut carry = 0;
        for i in (0..slots).rev() {
            let n = i as u32 + 2;
            let s = carry + tmp[i];
            let q = s / n;
            let r = s % n;
            tmp[i] = r * 10;
            carry = q;
        }
        digits[stop] = carry as u8;
    }
    digits
}

Back to π

After a slight detour, the study of π is resumed.