Alexander Bass
This note is one of many taken during 2026

Newton’s Method Pi

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

limngn=πwhereg0=3gn+1=gntan(gn) \begin{aligned} \lim_{n\to\infty} {g_n} &= \pi \\ &\text{where}\\ g_0 &= 3\\ g_{n+1} &= g_n - \tan(g_n) \end{aligned}

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

View Code
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 gn+1=gntan(gn)g_{n+1} = g_n - \tan(g_n) depends highly the initial value g0g_0. This makes some intuitive sense because sin(x)\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)=xtan(x)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 xx to the next value f(x)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 gn+1=gntan(gn)g_{n+1} = g_n - \tan(g_n) has a value which approaches the multiple of π nearest to g0g_0.

Specifically, let’s look at x=0π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 gn=gn+1g_n = g_{n+1}, or x=(xtan(x))x=-(x-tan(x)).

2x=tan(x) 2x = \tan(x)

Or

g(x)=tan(x)2x g(x)=\tan(x)-2x

Where g(x)=0g(x) = 0

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

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

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

g(x)=tan(x)2xg(x)=sec(x)tan(x)2 \begin{aligned} g(x) &= \tan(x) - 2x\\ g'(x) &= \sec(x)\tan(x) - 2 \end{aligned} vn+1=vntan(vn)2xsec(vn)tan(vn)2 v_{n+1} = v_n - \frac{\tan(v_n)-2x}{\sec(v_n)\tan(v_n)-2}

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

View Code
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π0\pi this well-behaved region is described as

x<S |x| < S

And in general is:

xnπ<S |x - n\pi| < S

Where:

nZS1.16556118 \begin{aligned} n &\in \mathbb{Z}\\ S &\approx1.16556118 \end{aligned}

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

Well, there’s the solution g0=Sg_{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 gng_n such that gn+1=gn+πg_{n+1}=g_{n}+\pi. A gng_n of this description will never converge.

xtan(x)+π=x x-\tan(x)+\pi = x tan(x)=π \tan(x) = \pi x=arctan(π) x = \arctan(-\pi)

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

arctan(±qπ) \arctan(\pm q\pi)

will diverge. For example:

g0=arctan(2π)g1=arctan(2π)2π \begin{aligned} g_0 &= \arctan(2\pi)\\ g_1 &= \arctan(2\pi)-2\pi \end{aligned}

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 π/2-\pi/2 (and all other midpoints between adjacent well-behaved areas)

As seen in the graph, the convergent values around π/2-\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).

View Desmos Code

Each line represents a distinct Desmos expression. Copy and paste them into Desmos to import it.

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

In the view, set the x-axis range to

π/20.1xπ/2+0.1 -\pi/2 - 0.1 \leq x \leq -\pi/2 +0.1

and the y-axis range to

100y100 -100 \leq y \leq 100

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

This note is one of many taken during 2026