Python Lab 1 part D

Newton's Method (NM)

The gradient descent method is the simplest version of a first order optimization method (a method which only uses first derivatives). Newton's method (NM) is a commonly used second order optimization method (a method which utilizes first and second derivatives). In this assignment you will implement NM, test its performance, and learn about its pitfalls.

1. Make a copy of the script you wrote for lab 1 part c, lj_1d.py, and rename it nm.py.
   
   
2. In NM, we will need to calculate second derivatives. Create a functions called second_deriv that calculates and returns the second derivative, V''(r_n), of the L-J potential.
   
   
3. Next change the name of the function gradient_descent to newtons_method and change your implementation of gradient descent in nm.py so that it implements NM. Here is how the positions should change with step number, n:
   r_n+1= r_n + F(r_n)/V''(r_n)
   where V''(r_n) is the second derivative of the potential.
   
   
4. To test your implementation, start your optimization at r=0.95. Once your code finds the analytical solution, you are done with this step.
   
   
5. Before moving on to the next step, think about the following question: How does NM performance compare to gradient descent?
   
   
6. Add a feature to the newtons_method function that exits out of your loop when you reach 100 optimization steps. There are multiple ways to do this, but for a new not introduced in prior python labs, google "python break statement".
   
   
7. Now, see if your optimizer finds the correct solution for the starting points r=1.2 and r=1.3. Note: One of these starting points should fail.
   
   
8. To understand why this is, create a few plots which test how the starting point of your optimization, r₀, affect the results your find with NM:
   
   
   1. Plot how the initial NM step size and direction (i.e. F(r_n)/V''(r_n)) changes with the initial position, r₀, from r = 0.9 to 1.5 with increments of 0.01.
      
      
   2. Plot how V''(r₀) changes with the initial position, r₀, from r₀ = 0.9 to 1.5 with increments of 0.01. It is particularly important to see when the curvature is negative or positive. To be able to clear this in your plot try setting the y axes range by adding the following line:
      
      ylim([-10.,10])
      
      You could also plot where the curvature is zero by plotting a horizontal line at zero curvature:
      
      plot([0.9,1.5],[0,0])
      
      
   3. Test whether your NM optimizer is successful for various starting positions, r₀ = 0.9 to r₀ = 1.5 with increments of 0.01. Plot whether the optimizer was successful (1) or not (0). Here are a few steps to make this happen:
      
      
      1. Make a function called test_nm. In this function create a for loop that goes through all your starting positions (r₀ = 0.9 to r₀ = 1.5 with increments of 0.01).
         
         
      2. Within the for loop you will need to optimize the L-J potential for various starting postions. Call newtons_method to do this in your for loop.
         
         
      3. You will also need to create an array indicating whether the optimizer was successful or not. Initialize an array prior to the for loop called success. Add a conditional statement within the for loop that ask whether the optimizer was successful. Depending on you result add a 0 or 1 to the array success. Remember when selecting the criteria in your conditional statement that in numerical calculations we get close to a critical point but do not find the exact value.
   
   
9. Looking at the plots you just made, try to answer the following:
   
   
   1. What happens to the initial NM steps when the optimizer first starts to fail moving from r₀=0.9 to larger values of r for the starting point?
      
      
   2. What happens when r₀ is greater than 1.25 (this value is approximate)?
   
   
   
10. What could we do to fix the problems in NM encountered in this lab so that we can find a local minimum at any starting point? If you have a good idea, try it and see how it performs!

Questions for 3C

Post all answers to a word document like Google Docs or Microsoft Word.

1. What is the path to your nm.py script?
2. Write down your answers to the following questions mentioned in step 8 above:
   1. What happens to the NM initial step when the optimizer first starts to fail and r is less than 1.25? Did it converge? Did the NM optimizer produce a local minimum? If it's not functioning as you expected, why did it turn out this way? (Optionally you can look at the trajectory of r, is it smooth? )
   2. What happens when r is greater than 1.25? Describe where the NM trajectory converges to, and whether or not your NM optimizer produced a local minimum.
   Include the plots you created in this lab to support your answers to these questions.
3. What is a potential solutions for fixing the two problems above? You can come up with this solution on your own or search the literature. Either way, let me know what the solutions is, how you came up with it, and how it fixes the problems.