# -*- coding: utf-8 -*-
"""
Created on Thu Jan 25 20:08:15 2022
@author: neil
"""

# program to illustrate iteration, slightly more robust than the previous
# this version has a convergence criteria, and also adapts the step size
# note that your starting point determines which zero you find!
# We iterate to find an x, where log(x)/2 == cos(x)
# Note this program uses our own function eqn(x), that returns cos(x)-log(x)/2

# (try using a seed like '10', which will illustrate that even this simple 
#  program can have major logic BUGS)

import numpy as np
import matplotlib.pyplot as plt

# function definition for the equation you are trying to find the zero
def eqn(x):                       # You can put any function here
    return np.cos(x)-np.log(x)/2    # return cos(x) - log(x)

# you need to select a starting point (this will determine which zero you find)
seedx = 4                   # starting point
# convergence criteria, how many sig digits?
conv = 0.00000001
dx=seedx/1000               # starting step size towards (hopefully) convergence
x = np.linspace(.2,7,100)   # set the 'grid' size for plotting (only)

# first plot what cos(x)-log(x) looks like (we could actually just look at the plot to see where eqq==0)
fig = plt.figure()
ax  = fig.add_subplot(111)
ax.plot(x,eqn(x),'b')       # plot residual of cos(x)-log(x), calls eqq()
ax.plot(x,np.cos(x),'g')    # plot log and cos functions
ax.plot(x,np.log(x)/2,'r')
ax.set_title('Log/2 and cos of x, and value of cos(x)- log(x)/2');
ax.set_xlabel('variable x')
ax.set_ylabel('values of functions')
ax.plot([0,7],[0,0],':k')   # draw a zero line on plot

# now start the iteration to get the program to find the zero of eqq()
# we are trying to reduce the residual from our starting point
rold = eqn(seedx)           # error or residual from initial guess
seedx += dx                 # start heading to right
rnew = eqn(seedx)           # residual from slightly different guess

while (np.abs(dx) > np.abs(conv)):
    if (np.abs(rold) < np.abs(rnew)):  # if our guess made it worse, change direction                           
        dx = -dx/2                     # make dx oposite sign, and take smaller steps
    seedx += dx                        # most of the direction changes are right around the zero
    rold = rnew                        # we have made a step, so our old position becomes rold
    rnew = eqn(seedx)                  # we keep going in dx steps until we cross zero and the sign
	                                   # changes, we then turn around and take 1/2 size steps, and repeat


ax.plot([seedx,seedx],[-1,1],':k')  # plot a vertical line at the zero point of the function
print("Numerical solution {}".format(seedx))
print("Approx Analytic solution to first root- 1.401")
plt.show()