# -*- coding: utf-8 -*-
"""
Created on Thu Jan 25 20:08:15 2021
Mainy problems require finding the minima or maxima of some function
or some other criteria.  A typical way of numerically attacking this 
is to try values in such a way that we get close to that criteria. We repeatedly 
chose values, if we are getting closer, we chose our next value in the same 
direction, if we are getting further away, we change direction.
In other words we iterate towards the solution.
@author: neil
"""

# program to illustrat iteration
# We iterate to find an x, where log(x)/2 == cos(x)
# a real program would at least have checks to make sure it doesn't
# get into an infinite loop
# Note this program uses our own function eqq(x), that returns cos(x)-log(x)/2
# the program is in 2 parts, the first is just plotting, the iteration is in the 2nd part

import numpy as np
import matplotlib.pyplot as plt

def eqq(x):                          # our function
    return np.cos(x)-np.log(x)*0.5   # return cos(x) - log(x)

seedx = 5                  # starting point
dx=seedx/10000            # (non adaptive)  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,eqq(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)*.5,'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 = eqq(seedx)           # error from initial guess
                            # take a small step and see if the error increases or decreases
rnew = eqq(seedx+dx)        # residual from slightly different guess
if np.abs(rnew)>np.abs(rold):
    dx = -1*dx              # if rnew is not smaller, then change the direction (change sign of dx)
seedx = seedx + dx          # dx should now be going in the right direction
rnew = eqq(seedx)

while np.abs(rnew)<np.abs(rold): # while the error is getting smaller
    rold = rnew
    seedx += dx
    rnew = eqq(seedx)

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 solutions  1.401, 5.783, 6.616")
plt.show()