10.4: Newton’s Method

We'll need the following in this section:

using Pkg; Pkg.add("Plots");

using Plots

Solving an equation is a very important part of mathematics and other scientific fields. However, finding general solutions to equations is not very easy. Consider a cubic equation like

\[15x^3 −143x^2 +226x+280=0 \label{cubic1} \] 

In the spirit of the quadratic formula, there is a cubic formula. Much of the wikipedia page spends time solving the cubic with all possibilities. In short, it’s not very easy.  In lieu of using such a formula, a more robust approach is to solve it numerically.  

This cubic equation on the left side of \ref{cubic1} actually factors, but finding those factors is quite difficult to do in general. A plot of the function is given with

f(x) = 15x^3-143x^2 +226x+280
plot([f,x->0],-2.5,10,label=["f(x)" "x-axis"])

The three intersection points between the red curve (\(x\)-axis) and the blue line (the function,\( y = f (x)\) ) are the three roots.

Newton’s method starts with a “guess” at the root and then refines it. Let \(x_0\) be the guess, then

We then look at using Newton’s method to solve this.  The following returns another value for the root.

\[x_1 = x_0 - \frac{f(x_0)}{f'(x_0)} \label{newton} \]

When Newton's method works (and often does well), the value of \(x_1\) is closer to the true root than \(x_0\).  We have the function \(f(x)\) defined in julia above, and can define the derivative with

df(x)=45x^2-286x+226

Waiting for kernel...

\[ df(x) = 45x^2-286x+226 \label{diff:cubic1} \]

We can then do a step of Newton's method with

x0=0
x1 = x0 - f(x0)/df(x0)

We can then do another step of Newton's method by using x1 as the input value and 

x2 = x1 - f(x1)/df(x1)

and note that it appears that the points are getting closer to the root (see the plot above) which is between 0 and -1.  A few more iterations of this should get a very accuracy value of the root. 

Exercise

Perform two more steps of Newton's method.  That is find x3 and x4.

# insert your code here.

A Newton's Method Function

It should appear natural to try to write a function to perform Newton's method for some number of steps.  We will walk through the thought process in developing the function.  We start with the template for this function.  The things that are needed to find a root using Newton's method is 1) the function 2) the derivative and 3) the initial guess.  Therefore we will start with

function newton(f::Function,df::Function,x0::Number) 

end

Here are considerations while building the function:

• You will need to do the two steps above many times so you will need a loop. Since you don’t know how many times you will need to run through the loop use a while loop and your condition should be that the two steps x0 and x1 are apart from each other.

• Checking if two floating-point numbers are equal are generally not a good idea, because they have to be equal to all bits, so instead we will run the while loop while the difference in the numbers x0 and x1 are larger than some default (like 10−6).

Here’s more a frame of the function:

function newton(f::Function, df::Function, x0::Number)
  x1 = x0 - f(x0)/df(x0)
  while abs(x1-x0)>1e-6 # while the two numbers are larger than 10^(-6)
    x0 = x1
    x1 = x0 - f(x0)/df(x0) 
  end
  x1
end

Waiting for kernel...

Using this we can now call this function as

newton(f,df,0)

Just to simplify this, we will define dx as -f(x0)/df(x0), which is the distance between successive x values, so we’ll use this to determine when to stop the loop:

function newton(f::Function, df::Function, x0::Real) 
  local dx=f(x0)/df(x0)
  while abs(dx)>1e-6
    x0 = x0-dx
    dx = f(x0)/df(x0)
  end
  x0
end

And just to ensure that the function still returns the same result, 

newton(f,df,0)

 

Using Automatic Differentiation

If you have used Computational Algebra Systems like Maple or Mathematica, you know that computers have the ability to differentiate. Julia does not have this capability built-in, although it has some capability of doing some of the feature set of these programs. There is a system called automatic differentiation that will compute the exact derivative to a function at a given point. That is, if you have a function f(x) and a number a, it will give you f′(a). The package is called ForwardDiff and you may need to add it and then

using Pkg; Pkg.add("ForwardDiff");
using ForwardDiff

For example, if you define: 

g(x) =x^2

g (generic function with 1 method)

ForwardDiff.derivative(g,3), ForwardDiff.derivative(g,-1) 

Newton’s Method with Automatic Differentiation

Using the ForwardDiff package, we can simplify Newton's method and not require that we need to find the derivative in calling it.  The following function 

function newton(f::Function, x0::Real)
  local dx = f(x0)/ForwardDiff.derivative(f,x0) 
  while abs(dx)>1e-6
    x0 = x0 - dx
    dx = f(x0)/ForwardDiff.derivative(f,x0)
  end
  x0
end

Waiting for kernel...