Exercise "Root Finding"

1. (6 points) Newton's method with numerical Jacobian and back-tracking linesearch

   1. Implement the Newton's method with simple backtracking linesearch algorithm where the derivatives in the Jacobian matrix are calculated numerically using finite differences.

      The interface to the function could be something like

      public static vector newton(
      	Func<vector,vector> f,
      	vector x,
      	double epsilon=1e-3,
      	double dx=1e-7
      ){
      	vector root;
      	/* calculate and */
      	return root;
      }
      

      where the parameters are:
      • function f: takes the input vector x and returns the vector f(x);
      • vector x: the starting point;
      • double epsilon: the accuracy goal: on exit the condition ‖f(x)‖<ε should be satisfied.
      • dx is the value of finite difference to be used in numerical evaluation of the Jacobian. It should be (much) smaller than the typical interval where the function changes appreciably but probably not smaller than square root of machine epsilon.

        One should stop one's iterations if either the condition ‖f(x)‖<ε is satisfied or if the step-size becomes smaller than the dx parameter.

        You should use your own routines for solving linear systems.

      • Find the extremum(s) of the Rosenbrock's valley function,

        f(x,y) = (1-x)2+100(y-x2)2
        

        by searching for the roots of its gradient.

2. (3 points) Bound states of hydrogen atom with shooting method for boundary value problems

   Introduction

   The s-wave radial Schrödinger equation for the Hydrogen atom reads (in units "Bohr radius" and "Hartree"),

   -(1/2)f'' -(1/r)f = εf ,
   

   where f(r) is the radial wave-function, ε is the energy, and primes denote the derivative over r.

   The bound s-state wave-function satisfies this equation and the two boundary conditions,

   f(r → 0) = r-r², (prove this)
   
   f(r → ∞) = 0 .
   

   These two boundary conditions can only be satisfied for certain discrete (negative) values of the energy.

   Since one cannot integrate numerically to ∞ one substitutes ∞ with a reasonably large number, r_max, such that it is much larger than the typical size of the hydrogen atom but still managable for the numerical inregrator (say, r_max = 10 Bohr radii),

   f(rmax)=0 .
   

   Let F_ε(r) be the solution (to be found numericall via gsl_odeiv) to our differential equation with energy ε and initial condition F_ε(r → 0)=r-r². Generally, for a random negative ε, this solution will not satisfy the boundary condition at r_max. It will only be satisfied when ε is equal one of the bound state energies of the system.

   Now define an auxiliary function

   M(ε) ≡ Fε(rmax) .
   

   The shooting method is then equivalent to finding the (negative) root of the equation

   M(ε) = 0 .
   

   Problems

   1. Find the lowest root, ε₀, of the equation M(ε) = 0 for, say, r_max=8. Plot the resulting function and compare with the exact result (which is ε₀=-1/2, f₀(r)=re^-r – check this by inserting ε₀ and f₀(r) into the Schredinger equation above).

   2. (Optional) Investigate the convergence of the solution towards the exact result as function of r_max.