Newton’s method

Newton’s method is technique to approximate solutions to equations that are difficult to solve numerically.

A simple applied problem

Here’s a simple applied problem, that doesn’t necessarily have a simple solution.

Statement

NoteMax area problem

We’d like to inscribe a rectangle with one vertex at the origin, opposite vertex on the graph of \(y=\sin(x)\) and base on the interval \([0,\pi]\). What rectangle meeting these conditions has maximal possible area?

The basic set up, which is not particularly hard, is shown in Figure 1. It’s pretty easy to see that we need to find the maximum value of \[ f(x) = x\sin(x) \] over the interval \([0,\pi]\). The image indicates that the maximum occurs a little bigger than \(x\approx2\), maybe? Can we find that exactly, though? If not, can we find at least find a more precise estimate?

viewof W = Inputs.range([0,3.141], 
    {value: 1.0, step:0.001, label: tex`W:`, format: d3.format("0.3f")}
)

plotOptions = ({
    width: Math.min(width, 600),
    height: 260,
    marginLeft: 45,
    marginRight: 20,
    x: {domain: [-0.1, 3.2], axis: null},
    y: {axis: null}
})

xAxis = () => Plot.axisX({
    y: 0,
    ticks: [Math.PI / 4, Math.PI / 2, 3 * Math.PI / 4, Math.PI],
    tickFormat: x => x === Math.PI / 2 ? "π/2" : x === Math.PI ? "π" : ""
})

curvePoints = d3.range(641).map(i => -0.1 + 3.3 * i / 640)

Plot.plot({
    ...plotOptions,
    y: {domain: [-1.2, 1.2], axis: null},
    marks: [
        Plot.rect([{x1: 0, x2: W, y1: Math.min(0, Math.sin(W)), y2: Math.max(0, Math.sin(W))}], {
            x1: "x1", x2: "x2", y1: "y1", y2: "y2",
            fill: "steelblue", fillOpacity: 0.2, stroke: "steelblue", clip: true
        }),
        Plot.ruleX([0]),
        Plot.ruleY([0]),
        xAxis(),
        Plot.axisY({x: 0, ticks: [-1, 1]}),
        Plot.line(curvePoints.map(x => [x, Math.sin(x)]), {strokeWidth: 2})
    ]
})

Plot.plot({
    ...plotOptions,
    y: {domain: [-1.2, 2.2], axis: null},
    marks: [
        Plot.ruleX([0]),
        Plot.ruleY([0]),
        xAxis(),
        Plot.axisY({x: 0, ticks: [-1, 2]}),
        Plot.line(curvePoints.map(x => [x, x * Math.sin(x)]), {strokeWidth: 1}),
        Plot.line(d3.range(401).map(i => {
            const x = W * i / 400;
            return [x, x * Math.sin(x)];
        }), {strokeWidth: 3, clip: true})
    ]
})

Figure 1: The graphs of \(y=\sin(x)\) and \(y=x\sin(x)\), with some decoration.

Aside on accuracy vs precision

Accuracy and precision are two related but distinct concepts in computational mathematics.

• An estimate is accurate if its actual error is small relative to the needs of a problem.
• An estimate is precise if its implied uncertainty is small.

As an example,

• \(\pi \approx 3.1\) is an accurate estimate of \(\pi\) with low precision,
• \(\pi \approx 3.4567891011\) is a precise estimate of \(\pi\) that’s not very accurate.

Generally, both are desirable but we should clear about which is which.

Start of a solution

It’s pretty easy to see that, as a function of \(x\), the area of the rectangle is \[ A(x) = x\sin(x). \] To use calculus to find the maximum possible area, we try to solve \(A'(x)=0\). That is \[ A'(x) = \sin(x) + x\cos(x) \stackrel{?}{=} 0. \]

You can have a go at that, if you’d like!

Solving the equation on a computer

These days, I guess the first computational tool you might think to try could be a chatbot. Here’s a transcription of Chat GPT’s response to my query on this question:

ChatGPT’s solution

1. You said:

2. ChatGPT said:

   We want to solve

   $$\sin x+x\cos x=0.$$

   Step 1: Rearrange the equation

   Move one term to the other side:

   $$\sin x=-x\cos x.$$

   If $\cos x\neq 0$, divide both sides by $\cos x$:

   $$\tan x=-x.$$

   Therefore, we need to solve

   $$\boxed{\tan x+x=0}.$$

   Step 2: Find the solutions

   Notice that $x=0$ is a solution because

   $$\sin 0+0\cos 0=0.$$

   To determine whether there are other solutions, observe that the original equation is equivalent to

   $$\sin x+x\cos x=0.$$

   For $x\neq 0$, this can be written as

   $$\frac{\tan x}{x}=-1.$$

   The equation is transcendental, so its solutions generally cannot be expressed in closed form using elementary functions.

   Also, the original equation is odd, meaning that if $x$ is a solution, then $-x$ is a solution.

   The nonzero solutions can be found numerically. They satisfy

   $$\tan x=-x.$$

   For example, the first positive nonzero solution is approximately

   $$x\approx 2.029.$$

   By symmetry, another solution is $x\approx -2.029$.

   Final answer

   The solutions include

   $$\boxed{x=0,\qquad x\approx\pm 2.029,\ \pm 4.913,\ \pm 7.978,\ldots}$$

   These are numerical solutions; there are infinitely many real roots.

3. You said:

4. ChatGPT said:

   Good question! Unfortunately, there isn't a known exact answer in the form of simple numbers like $\frac12$, $\sqrt2$, or $\pi$ for the nonzero solutions.

   We have

   $$\sin x+x\cos x=0$$

   Dividing by $\cos x$, when $\cos x\ne0$, gives

   $$\tan x=-x.$$

   Why can't we solve it exactly?

   For example, if we try $x=\frac{\pi}{2}$, then $\cos x=0$, but $\sin x=1$, so it doesn't work.

   If we try $x=\sqrt2$, we'd need

   $$\tan(\sqrt2)=-\sqrt2,$$

   which isn't true.

   The difficulty is that $x$ appears both inside a trigonometric function and outside it. This makes the equation transcendental, so standard algebraic techniques generally can't isolate $x$.

   What is the exact answer, then?

   We can write the exact solution set implicitly as

   $$\boxed{x=0\quad\text{or}\quad \tan x=-x.}$$

   The nonzero roots must be approximated numerically, for example,

   $$x\approx\pm2.029,\quad \pm4.913,\quad \pm7.978,\ldots$$

   One important distinction: not having a simple exact form doesn't prove that a root cannot be expressed using some more advanced mathematical function. It means there is no standard elementary closed-form solution.

Using SageMath

It’s not hard and worth understanding how to do this sort of thing with a real computer algebra system, though. There are lots of tools for this, including Mathematica, Matlab, Python, Sympy, and SageMath. I’ve used all of these over the years.

One nice aspect of SageMath is that it can easily be embedded into a webpage. Here’s how to use SageMath to approach this problem, both algebraically, and numerically:

It turns out that SageMath gives us 16 digits of precision. Generally, I would expect almost all of those to be accurate.

Newton’s algorithm

Newton’s method is an algorithm to find numerical approximations to roots of functions; it is the theoretical foundation on which numerical tools like Sage’s find_root works. Given an initial guesss \(x_1\), Newton’s method improves this guess by applying the function \[N(x) = x - \frac{f(x)}{f'(x)}.\] This produces \(x_2 = N(x_1)\). We then plug that back in to get \(x_3\) and continue. More generally, we produce a sequence \((x_k)\) via \(x_k=N(x_{k-1})\).

Example

Let \(f(x) = x^3-x-1\). It’s evident from a graph that there’s one root a little larger than 1.

{
    const f = x => x**3 - x - 1;
    const pts = d3.range(-2.1,2.8,0.01).map(x => [x,f(x)])
    return Plot.plot({
        marks: [
            Plot.line(pts, {strokeWidth: 3}),
            Plot.ruleX([0]), Plot.ruleY([0]),
            Plot.axisX({y:0, ticks: [-2,-1,1,2]}), 
            Plot.axisY({x:0, ticks: [-8,-2,4,8,12,16]})
        ]
    })
}

Newton’s method works by riding the tangent line from an initial guess. If we note that \(x_1=2\) is pretty close to the root, we compute \(x_2 = N(x_2)\), where \[N(x) = x - \frac{f(x)}{f'(x)} = x - \frac{x^3-x-1}{3x^2-1}.\] Thus, \(x_2 = N(2) = 2-5/11 \approx 1.54545\). Geometrically, this point is obtained by riding the tangent line to the \(x\)-axis:

{
    const f = x => x**3 - x - 1;
    const fp= x => 3*x**2 - 1;
    const x1 = 2;
    const y1 = f(x1);
    const m1 = fp(x1);
    const l1 = x => m1*(x-x1)+y1;
    const x2 = x1 - f(x1)/fp(x1);
    const pts = d3.range(-2.1,2.8,0.01).map(x => [x,f(x)]);
    return Plot.plot({
        y: {domain: [-4,17]},
        marks: [
            Plot.line(pts, {strokeWidth: 3}),
            Plot.line([[0,l1(0)], [3,l1(3)]], {strokeWidth: 1}),
            // Plot.line([[0,l2(0)], [3,l2(3)]], {strokeWidth: 1}),
            Plot.dot([[x1,y1]], {fill: "green", r: 6}),
            Plot.dot([[x2,0]], {fill: "red", r: 4}),
            Plot.ruleX([0]), Plot.ruleY([0]),
            Plot.axisX({y:0, ticks: [-2,-1,1,2]}), 
            Plot.axisY({x:0, ticks: [-2,4,8,12,16]})
        ]
    })
}

If we do that again, we end up even closer to the root:

{
    const f = x => x**3 - x - 1;
    const fp= x => 3*x**2 - 1;
    const x1 = 2;
    const y1 = f(x1);
    const m1 = fp(x1);
    const l1 = x => m1*(x-x1)+y1;
    const x2 = x1 - f(x1)/fp(x1);
    const y2 = f(x2);
    const m2 = fp(x2);
    const l2 = x => m2*(x-x2)+y2;
    const x3 = x2 - f(x2)/fp(x2);
    const pts = d3.range(-2.1,2.8,0.01).map(x => [x,f(x)]);
    return Plot.plot({
        y: {domain: [-4,17]},
        marks: [
            Plot.line(pts, {strokeWidth: 3}),
            Plot.line([[0,l1(0)], [3,l1(3)]], {strokeWidth: 1}),
            Plot.line([[0,l2(0)], [3,l2(3)]], {strokeWidth: 1}),
            Plot.dot([[x1,y1], [x2,f(x2)]], {fill: "green", r: 6}),
            Plot.dot([[x2,0], [x3,0]], {fill: "red", r: 4}),
            Plot.ruleX([0]), Plot.ruleY([0]),
            Plot.axisX({y:0, ticks: [-2,-1,1,2]}), 
            Plot.axisY({x:0, ticks: [-2,4,8,12,16]})
        ]
    })
}

That’s why we iterate!

Performing Newton’s method with Sage

Here’s how we can apply Newton’s method to the previous example using Sage:

Performing Newton’s method by hand

A classic problem is to use Newton’s method to approximate \(\sqrt{2}\) by hand. If we stick with fractions, we will generate good rational approximations to \(\sqrt{2}\). The result is equivalent to an algorithm known by the Babylonians.

To derive the technique using Newton’s method, define \(f\) by \[ f(x) = x^2 - 2. \] Clearly, the positive root of \(f\) is exactly \(\sqrt{2}\). Now, \[ N(x) = x - \frac{f(x)}{f'(x)} = x - \frac{x^2-2}{2x} = \frac{x}{2} + \frac{1}{x}. \] Now, let’s take two Newton steps from \(x_0=1\): \[ \begin{aligned} x_1 &= \frac{x_0}{2} + \frac{1}{x_0} = \frac{1}{2} + \frac{1}{1} = \frac{3}{2} \\ x_2 &= \frac{x_1}{2} + \frac{1}{x_1} = \frac{3/2}{2} + \frac{1}{3/2} = \frac{3}{4} + \frac{2}{3} = \frac{17}{12}. \\ \end{aligned} \] Note that \[ \left(\frac{17}{12}\right)^2 = \frac{289}{144} = \frac{288}{144} + \frac{1}{144} = 2 + \frac{1}{144} \] is really close to two! So, \(17/12\) must be really close to \(\sqrt{2}\).

Exercises

1. Use a numerical tool to solve the following equations. Be sure to find all solutions
   1. \(x^5-x-1 = 0\)
   2. \(x^5-2x-1 = 0\)
   3. \(\sin(3x) = x/2\)
2. For each of the following functons, take three Newton steps from the given initial point
   1. \(f(x) = x^2 - 2\) from \(x_1 = 2.0\)
   2. \(f(x) = \sin(x)\) from \(x_1 = 3.0\)
   3. \(f(x) = x^5-x-1\) from \(x_1 = 1.0\)
3. The equation \(\sin(x)=x/9\) has 7 solutions. We want to find an approximation to the largest solution. Use a graph to find a good initial approximation for your \(x_1\) and apply Newton’s method from that point obtain the approximation.