#!/usr/bin/env python
# coding: utf-8

# # Finding roots of equations - Newton-Raphson method and secant method
# 
# In this lesson, we will find the roots of an equation using two iterative methods. By the end of the lesson, you will be able to:
# 1.	describe what a root of an equation is.
# 2.	approximate the roots of an equation using the Newton-Raphson method.
# 3.	approximate the roots of an equation using the secant method.
# 
# ## Roots of an equation
# What is a root? Consider a function $f(x)=x^2-5x+6$. A root of the equation $f(x)$ is a value $x_0$ such that $f(x_0)=0$. In this case, $f(x)$ has two roots: $x=2$ and $x=3$.
# $$f(2)=2^2-5 \times 2 + 6 = 0$$ 
#  
# If we plot the function   we can see where the graph crosses the x-axis. The points where it crosses the x-axis are the roots of $f(x)$. Don't worry about the code used to plot the function below. We will learn how to plot graphs later. Just run the cell.

# In[ ]:


import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(1, 3.5, 100)
plt.figure()
plt.subplots()
plt.plot(x, x**2-5*x+6)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(b=True)
plt.show()


# For the above function, the roots $x=2$ and $x=3$ can be found analytically. However for many situations, the roots of an equation cannot be determined analytically so numerical methods for estimating them have been developed. Newton's method is one such approach. Newton's method is frequently used in a variety engineering applications. These would include optimising processes to find the best operating parameters or determining the material parameters that describe the mechanical behaviour of complex material.
# 
# ## Newton-Raphson method
# The Newton-Raphson method or Newton's method is a technique for solving equations numerically. Specifically, it approximates the roots of equations by making a guess at the root and then repeatedly updating the guess to arrive at a satisfactory approximation of the root.
# 
# Consider a function $f(x)$ with a root at **o** as indicated in the diagram. Newton's method takes the following steps to estimate **o**:
# Guess an initial value for the root $x_0$ .
# 1.	The slope at $f(x_0)$ is calculated. Where the slope crosses the x-axis is an updated guess for the root $x_1$.
# 2.	The slope at $f(x_1)$ is calculated. Where the slope crosses the x-axis is an updated guess for the root $x_2$ .
# 3.	This procedure is repeated until the root of $f(x)$ is determined within a specified tolerance.
# 
# Note in the diagram how $x_1$ is closer to **o** than $x_0$ and that $x_2$ is closer than $x_1$.
# 
# <img src="./img/newtons_method.svg.png" alt="Newton's method plotted" align='center' style="width: 600px;"/>
# 
# It can be shown that $x_1$ is given by
# $$x_1=x_0-\frac{f\left(x_0\right)}{f^\prime\left(x_0\right)}$$
# 
# where $f^\prime\left(x_0\right)$ is the slope or first derivative of $f(x)$ at $x_0$.
# 
# This can be generalised to
# $$x_{n+1}=x_n-\frac{f\left(x_n\right)}{f^\prime\left(x_n\right)}$$
# We are going to create a basic Python script that uses Newton's method to find the root of an equation. The first script will be just enough to find the root. We will then develop the basic script to output results and messages.

# In[ ]:


# Function that we are finding the root of
def fx(x):
    
    y = x**2 - 5*x + 6
    
    return y

# Derivative of the function that we are finding the root of
def dfdx(x):
    
    y = 2*x - 5

    return y

# This script calculates the roots of an equation
x0 = 1.0 # This is the initial guess
max_iter = 20 # % This is the maximum number of iterations allowed
tolerance = 1e-4 #% Tolerance for which f(x) must be less than to determine root

num_decimal_places = 4

for i in range(1,max_iter + 1):

    if abs(fx(x0)) < tolerance:
        print(f'\nThe root of the equation is {round(x0,num_decimal_places)}.\n')
        break
    
    x1 = x0 - fx(x0) / dfdx(x0); # Newton's method equation
    x0 = x1; # Updating the guess


# Note that you must define two functions - one to define $f(x)$ and another to define the derivative of the function $f^{\prime}(x)$.
# 
# In the code above, change the value of tolerance to a bigger or smaller value. What is the effect on the result?
# 
# ## Exercises
# Modify the code above to incorporate the following features:
# 1.	Output the results of each iteration in the form of a table. The first column gives the iteration number, the second column gives $x_0$, and the third column gives  $f(x_0)$.
# 2.  Put the code for the Newton method in a function. The input arguments for the function are the initial guess, maximum iterations, and the tolerance. The output arguments are a variable indicating whether the root is found and the root of the equation.
# 3.  Call the function in 2 to find the root of the equation. Display appropriate messages to the shell indicating whether the root has been found and what it is.
# 3.	Use a `while` loop in place of a `for` loop.
# 
# ### Problem 1
# A cantilever beam is loaded with a distributed load. The deflection y of the centreline of the beam is a function of the position x:
# $$y=f(x)=\frac{w_oL}{3\pi^4EI}\left(48L^3cos\left(\frac{\pi}{2L}x\right)-48L^3+3\pi^3Lx^2-\pi^3x^3\right)$$
# where $L$ is the length of the cantilever, $E$ is the elastic modulus, $I$ is the second moment of area, and $W$ is the amplitude of the applied load.
# Let $L=4$ m, $E=70$ GPa, $I=52.9 \times 10^{-6}$ $\text{m}^4$, $w_o=25$ kN/m.
# Using Newton's method, find the positions $x$ to four decimal places when the deflection of the beam $y=5$ mm.
# Use an initial guess of $x=2$ m and a tolerance value of 0.0001.
# 
# ### Problem 2
# The variation in temperature of a copper rod when a current is applied is given by the following first order differential equation:
# $$\frac{dT}{dt}=\frac{I^2R^{\prime}_e-\pi Dh\left(T-T_{\alpha}\right)-\pi D \epsilon \sigma \left(T^4 - T^4_{sur}\right)}{\rho c \left(\pi D^2/4\right)}$$
#  
# where $I=10$ A, $R^{\prime}_e=0.4$ $\frac{\Omega}{\text{m}}$, $T_{sur}=300$ K, $T_{\alpha}=300$ K, $h=100$ $\frac{\text{W}}{\text{m}^2}$, $\epsilon=0.8$, $\sigma=5.67 \times 10^{-8}$ $\frac{\text{W}}{\text{m}^2\text{K}}$, $D=0.001$ m, $\rho=8800$ $\frac{\text{kg}}{\text{m}^3}$, $c=420$ $\frac{\text{J}}{\text{kgK}}$.
# 
# Find the steady state temperature values in kelvin by using the Newton method. Use an initial guess of 300 K and apply a tolerance value of 0.001.
# 
# ## Secant method
# The secant method is similar to Newton's method. Instead of the derivative to the function, an approximation to the derivative is provided. Two initial guesses are provided instead of one. The formula for updating the estimate to the root   is
# $$x_{n+1}=x_n-f\left(x_n\right)\frac{x_n-x_{n-1}}{f\left(x_n\right)-f\left(x_{n-1}\right)}$$
#  
# The following script finds a root of the equation $x^2-5x+6$ using the secant method. Copy the following code into a Python script file and run it.

# In[ ]:


def gx(x):
    
    y = x**2 - 5*x + 6
    
    return y

# This script calculates the roots of an equation using the secant method

x0 = 4 # This is the first guess
x1 = 4.5 # This is the second guess
max_iter = 20 # This is the maximum number of iterations allowed
tolerance = 1e-4 # Tolerance for which f(x) must be less than to determine root

for i in range(1,max_iter + 1):

    if abs(gx(x0)) < tolerance:
        print(f'\nThe root of the equation is {round(x0,4)}.\n')
        break
    elif abs(gx(x1)) < tolerance:
        print(f'\nThe root of the equation is {round(x1,4)}.\n')
        break

    x2 = x1 - gx(x1) * (x1 -x0) / (gx(x1) - gx(x0)); # Secant method equation to update guess
    x0 = x1;
    x1 = x2;


# ## Exercises
# Modify the code above to incorporate the following features:
# 5.	Output the results of each iteration in the form of a table. The first column gives the iteration number, the second column gives $x_0$, the third column gives  $x_1$, the fourth column gives $g(x_0)$, and the fifth column gives $g(x_1)$.
# 6.	If the root is not found within the maximum number of iterations, a message is output to the command window.
# 
# ## Problem 3
# The viscous or frictional loss in pipe flow expressed in terms of the friction factor $f$ is given by Colebrook’s equation:
# 
# $$\frac{1}{f^{0.5}}=-2\ln\left(\frac{\epsilon/D}{3.7}+\frac{2.51}{Ref^{0.5}}\right)$$
#  
# where $\epsilon=90e-6$ m is the roughness, $D=0.15$ m is the pipe diameter, $Re=3 \times 10^5$  is the Reynold's number.
# Find the friction factor $f$ by the secant method. Apply a tolerance value of 0.001 with initial guess values of 0.006 and 0.005.
# 
# ## Problem 4
# A soccer player is about to do a long throw into his striker. The footballer is 1.7 m tall and is supposed to hit the striker 30 m away. The equation that describes the motion of the football is derived from projectile motion and is given by
# 
# $$y=x \tan{\theta}-\frac{1}{2}\frac{x^2g}{v_o^2}\frac{1}{\cos^2{\theta}}+h$$
#  
# where $x=30$ m and $y=2$ m are the horizontal and vertical distances, respectively; $g$ is the acceleration due to gravity. $v_0=25$ m/s is the initial velocity of the ball as it leaves the footballer's hands, $h=1.7$ m is the height of the footballer and $\theta$ is the angle the ball makes with the horizontal just as it leaves the footballer's hands.
# 
# Find the angle $\theta$ that the footballer must launch the ball by the secant method to four decimal places. Apply initial guesses of 0.6 and 0.65 rad. Set the tolerance at 0.0001.

# In[ ]: