How to Find a Root of a Polynomial using Newton Raphson Method

Newton-Raphson Method of finding roots of a polynomial is a Householder Method of order 1.Householder's Methods are used to find roots for  functions of one real variable with continuous derivatives up to some order. Here is how this works for order 1.

• We set an approximate value for the root (x0). We take the value of the function (f) and its derivative(g) at this point.
• The derivative at this point is simply the tangent to the function at this point. Now, we find out where this tangent intersects the x-axis. We set this point as the new x0.
• We go on calculating the function and its derivative at each x0 we encounter for 'm' number of iterations. Finally. we reach a point where the function attains zero value and f(x0)/g(x0) is zero. This gives us root=the current value of x0.

Pros

• Simple Logic: Easy to Code
•  Quadratic Convergence:  The number of accurate digits double with each iteration.
• Works best when you set the number of iterations and work within an interval where a function and its derivative is well-behaved.

Cons:

• Requires knowledge of the derivative, which can slow down calculations.
• Requires assumptions in the proof of quadratic convergence to be true.
• May encounter an overshoot (deviation from the root in a neighborhood where the derivative is not well-behaved), or a stationary point (where the derivative is 0 and f/g cannot be calculated).
• Requires additional steps (successive over-relaxation) to maintain quadratic convergence rate for multiplicity other than unity.

Now for the FORTRAN program,

Note: I have added a subroutine to generate a polynomial function. However, you can simply insert a function and its derivative. I found this easier.

 !To find the  root of a polynomial using Newton Raphson Method
program newt_raph
    implicit none
    integer iter,n,m,i,n1
    real x0,root,sum1,sum0,g,f,eps
    real, allocatable :: a(:)

    write(*,*) "Give the degree of the polynomial" !To read the coefficients of the polynomial function
    read(*,*) n1
    write(*,*) "For a polynomial of degree, ",n1

    n=n1+1
    allocate (a(0:n))
    write(*,*) "Give the coefficients of the polynomial"
    do i=n-1,0,-1
         read (*,*) a(i)
      end do

      write (*,*) "The coefficients of the polynomial are"
      do i=n-1,0,-1
        write(*,*) a(i)
    end do

       write(*,*) "Give the approximate value of the root" !To set an initial guess value for the root
    read (*,*) x0
    write(*,*) "The approximate root is ",x0

    write (*,*) "Give the number of iterations" !For Newton Raphson Method
    read (*,*) m
    write (*,*) "Maximum ",m," iterations to be carried out"

    sum1=x0
    iter=0
    eps=1.0e-5

    do
          sum0=sum1 !To avoid division by zero
          call horner (a,n,sum1,f,g)
          if (abs(g)<eps) then
            write (*,*) "Derivative too small"
               stop
        end if
      
          sum1=sum1-(f/g) !To calculate the root
          iter=iter+1
        if ((abs(sum1-sum0))<0.00001) then
              root=sum1
              write (*,*) "The root by Newton Raphson Methos is ",root," ising ",iter," iterations"
              exit
          end if
      
          if (iter>m) then !To exit on not being able to find the root
            write (*,*) "Not able to converge with the given number of iterations"
            exit
        end if
      
      end do

end program newt_raph

    subroutine Horner (a,n,x,p,q) !To calculate the function using Horner's Method (optional)
        implicit none
        integer i,n
        real p,q,x,b(0:n),a(0:n)
        p=a(n-1) !For the function
        do i=n-2,0,-1
              p=p*x+a(i)
          end do
          do i=(n-2),0,-1 !For the derivative
               b(i)=a(i+1)*(i+1)
        end do
        q=b(n-2)
        do i=n-3,0,-1
              q=q*x+b(i)
          end do
    end subroutine Horner

Sample Output for Newton Raphson Method

Gnuplot depiction of the Roots of the Function