# Root finding (Iteration)

Below program use simple iteration to find the root of a function of the form of
\begin{equation}
f(x) = x - rx\left(1+\dfrac{x}{4}-x^3\right)
\tag{1}
\end{equation}

First, I use the following iterative equation.
\begin{equation}
x = rx\left(1+\dfrac{x}{4}-x^3\right)
\tag{2}
\end{equation}

I find out that $x$ doesn't settle to a single value most of the time. Futhermore, it became more diverge as the value of $r$ increase. Only when $r=0.8$ it slowly converge to $0$.


If I then use another iterative equation, it diverge super fast to infinity.
\begin{equation}
x = \dfrac{r}{1-r}+rx^2\left(\dfrac{1}{4}-x^2\right) \tag{3}
\end{equation}

## r=0.8

## r=1.2

I think the reason lies in the fact that this root finding method is very sansitive to its initial value and iterative formulation. As you can see from the graph of $f(x)$, as $r$ increase the slope around the bottom of the graph became steeper around $x=0$ makes the iteration less likely to converge to a single value.

## Code

!========================================================
!   Purpose: Find root
!
!   Methods: naive iteration
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!     10/15/20       MorrisH        Original Code
!========================================================
PROGRAM MAIN
    IMPLICIT NONE
    REAL(8) :: r, x
    INTEGER :: i

    open(66, file='20.dat')
    CALL RANDOM_SEED()
    CALL RANDOM_NUMBER(x)
    r = 2.0d0
    do i = 1, 100
        write(66, *) x
        x = r * x * (1.d0 + x / 4.d0 - x**3) 
    end do
END PROGRAM MAIN

# (Real) Root finding
Following program devide into four parts:
* main : Integrate all subroutines and analysis CPU\_TIME.
* BISECTION : Use bisection method to evaluate intergral.
* GOLDENBISECTION : Change bisection by golden-ratio-section to evaluate integral.
* NewtonRaphson : Use Newton-Raphson to evaluate the integral.
\begin{equation}
x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)} \tag{4}
\end{equation}

Testing function
\begin{equation}
f(x) = 2-e^{-x^2/8}\left(x^2-x+1\right)
\tag{5}
\end{equation}
As you can see from above lsit, there are four solution(root) to this function.
I use three seperate method to compare it's time and accuracy. Both three method are roughly the same accuracy(only different by around $10^{-12}$). But the Newton-Raphson takes less time and much less steps in general. Two slightly different bisection method takes roughly the same time and steps to finish. One interesting thing to notice is that Newton-Raphson method are very sensitive to it's initial value, slighly differnt can result into to
diverge the solution.


## Result


### root 1


### root 2

### root 3

### root 4

## Code

### Main

PROGRAM MAIN
    IMPLICIT NONE
    REAL(8) :: f, a, b, eps, x, start, finish
    INTEGER :: n
    EXTERNAL :: f, fd

    a = 3.1d0
    b = 4.d0
    eps = 1.e-12
    print *, "Finding root between", a, "and", b
    print *, "with the error be less than:", eps
    write(*, '(a)') repeat('-', 70)

    CALL CPU_TIME(start)
    CALL BISECTION(f, a, b, eps, x, n)
    CALL CPU_TIME(finish)
    print *, "Bisection method:"
    print *, "root:", x
    print *, "takes ", n, "step"
    print *, "takes ", finish-start, "sec"
    write(*, '(a)') repeat('-', 70)

    CALL CPU_TIME(start)
    CALL GoldenBISECTION(f, a, b, eps, x, n)
    CALL CPU_TIME(finish)
    print *, "Bisection method(alt):"
    print *, "root:", x
    print *, "takes ", n, "step"
    print *, "takes ", finish-start, "sec"
    write(*, '(a)') repeat('-', 70)

    CALL CPU_TIME(start)
    CALL NewtonRaphson(f, fd, a, eps, x, n)
    CALL CPU_TIME(finish)
    print *, "Newton-Raphson method:"
    print *, "root:", x
    print *, "takes ", n, "step"
    print *, "takes ", finish-start, "sec"
    write(*, '(a)') repeat('-', 70)
END PROGRAM MAIN

### Bisection

SUBROUTINE BISECTION(f, a, b, eps, x, n)
    !========================================================
    !   Purpose: Find root
    !
    !   Methods: Bisection method
    !       
    !       Date        Programer       Description of change
    !       ====        =========       =====================
    !     10/15/20       MorrisH        Original Code
    !--------------------------------------------------------
    !   input:
    !           func -- function to evaluate           
    !           a    -- xmin
    !           b    -- xmax
    !           eps  -- accuracy
    !   output:
    !           x    -- root
    !           n    -- iterative times
    !========================================================
    IMPLICIT NONE
    REAL(8), EXTERNAL :: f
    REAL(8) :: a, b, eps, x
    REAL(8) :: x1, x2, x3
    INTEGER :: i, n, nmax=1000

    ! check bisection condition
    if ( f(a) * f(b) > 0.d0 ) then
        print *, "There is no root between a and b"
        return
    end if

    x1 = a
    x3 = b
    do i = 1, nmax

        x2 = (x1 + x3) / 2.d0 
        if ( f(x1) * f(x2) < 0.d0 ) then
            x3 = x2
        else
            x1 = x2
        end if

        if ( abs(x1-x3) <= eps ) exit
    end do
    n = i
    x = (x1 + x3) / 2.d0

END SUBROUTINE BISECTION

### GoldenBisection

SUBROUTINE GoldenBISECTION(f, a, b, eps, x, n)
    !========================================================
    !   Purpose: Find root
    !
    !   Methods: Bisection method + divide by golden mean
    !       
    !       Date        Programer       Description of change
    !       ====        =========       =====================
    !     10/15/20       MorrisH        Original Code
    !--------------------------------------------------------
    !   input:
    !           func -- function to evaluate           
    !           a    -- xmin
    !           b    -- xmax
    !           eps  -- accuracy
    !   output:
    !           x    -- root
    !           n    -- iterative times
    !========================================================
    IMPLICIT NONE
    REAL(8), EXTERNAL :: f
    REAL(8) :: a, b, eps, x
    REAL(8) :: x1, x2, x3
    INTEGER :: i, n, nmax=1000

    ! check bisection condition
    if ( f(a) * f(b) > 0.d0 ) then
        print *, "There is no root between a and b"
        return
    end if

    x1 = a
    x3 = b
    do i = 1, nmax

        x2 = x1 + (x3 - x1) * (dsqrt(5.d0) - 1) / 2 
        if ( f(x1) * f(x2) < 0.d0 ) then
            x3 = x2
        else
            x1 = x2
        end if

        if ( abs(x1-x3) <= eps ) exit
    end do
    n = i
    x = x1 + (x3 - x1) * (dsqrt(5.d0) - 1) / 2 


END SUBROUTINE GoldenBISECTION

### NewtonRaphson

\subsection*{NewtonRaphson.f90}
\begin{minted}[frame=lines, linenos, fontsize=\large]
{fortran}
SUBROUTINE NewtonRaphson(f, fd, a, eps, x, n)
    IMPLICIT NONE
    REAL(8), EXTERNAL :: f, fd
    REAL(8) :: a, eps, x, old
    INTEGER :: i, n, nmax=1000

    x = a
    do i = 1, nmax
        old = x
        x = x - f(x) / fd(x)
        if ( abs(old - x) <= eps ) exit
    end do
    n = i
END SUBROUTINE NewtonRaphson
\end{minted}


\subsection*{func.f90}
\begin{minted}[frame=lines, linenos, fontsize=\large]
{fortran}
FUNCTION f(x)
    IMPLICIT NONE
    REAL(8) :: f, x
    f = 2.d0 - exp(-x * x / 8.d0) * (x * x - x + 1)
    ! f = x
END FUNCTION f


FUNCTION fd(x)
    IMPLICIT NONE
    REAL(8) :: fd, x
    fd = exp(-x * x / 8.d0) * (x**3 - x**2 - 7.d0 * x + &
                               4.d0) / 4.d0
    ! fd = 1.d0
END FUNCTION fd