Homework "Root Finding"

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

   Implement the Newton's method with simple backtracking line-search algorithm where the derivatives in the Jacobian matrix are calculated numerically using finite differences as described in the Book.

   Something like this,

   static vector newton(
   	Func<vector,vector>f /* the function to find the root of */
   	,vector start        /* the start point */
   	,double acc=1e-2     /* accuracy goal: on exit ‖f(x)‖ should be <acc */
   	,vector δx=null      /* optional δx-vector for calculation of jacobian */
   	){
   vector x=start.copy();
   vector fx=f(x),z,fz;
   do{ /* Newton's iterations */
   	if(fx.norm() < acc) break; /* job done */
   	matrix J=jacobian(f,x,fx,δx);
   	var QRofJ = givensQR(J);
   	vector Dx = QRofJ.solve(-fx); /* Newton's step */
   	double λ=1;
   	do{ /* linesearch */
   		z=x+λ*Dx;
   		fz=f(z);
   		if( fz.norm() < (1-λ/2)*fx.norm() ) break;
   		if( λ < λmin ) break;
   		λ/=2;
   		}while(true);
   	x=z; fx=fz;
   	}while(true);
   return x;
   }
   

   The routine returns a vector x that approximates the root of the equation f(x)=0 such that ‖f(x)‖<acc.

   The vector δx to be used in the finite-difference numerical evaluation of the Jacobian depends on the problem at hand and should be supplied by the user (who should know their problem better). If the user does not supply it, the routine can choose it as

   δxi = |xi|*2-26 ,
   

   or as (might work better sometimes),

   δxi = Max(|xi|,1)*2-26 ,
   

   where 2^-26 is the square root of the (double) machine epsilon.

   The Jacobian can be estimated numerically like this,

   public static matrix jacobian
   (Func<vector,vector> f,vector x,vector fx=null,vector dx=null){
   	if(dx == null) dx = x.map(xi => Abs(xi)*Pow(2,-26);
   	if(fx == null) fx = f(x);
   	matrix J=new matrix(x.size);
   	for(int j=0;j < x.size;j++){
   		x[j]+=dx[j];
   		vector df=f(x)-fx;
   		for(int i=0;i < x.size;i++) J[i,j]=df[i]/dx[j];
   		x[j]-=dx[j];
   		}
   	return J;
   }
   

   One should stop one's iterations if either the condition ‖f(x)‖<acc is satisfied (or if the step-size becomes smaller than the size of the δx parameter in your numerical gradient calculation).

   You should use your own routines for solving linear systems.

   Debug your root-finding routine using some simple one- and two-dimensional equations.

   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 (you should calculate the latter analytically).

   Find the minimum(s) of the Himmelblau's function,

   f(x,y) = (x2+y-11)2+(x+y2-7)2 ,
   

   by searching for the roots of its (analytic) 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 = Ef ,
   

   where f(r) is the (reduced) radial wave-function, E 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 the Coulomb potential diverges at zero one cannot specify the boundary condition at r=0 numerically. Instead one can use some finite radius r_min, that is much smaller than the Bohr radius, and specify the boundary condition at r_min,

   f(0)=0 → f(rmin)=rmin-rmin2 .
   

   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 manageable for the numerical integrator (say, r_max = 10 Bohr radii),

   f(∞)=0 → f(rmax)=0 .
   

   Let F_E(r) be the solution (to be found numerically using your own ODE solver) to our differential equation with energy E and the initial condition F_E(r_min)=r_min-r_min². Generally, for a random negative E, this solution will not satisfy the boundary condition at r_max. It will only be satisfied when E is equal one of the bound state energies of the system.

   Now define an auxiliary function

   M(E) ≡ FE(rmax) .
   

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

   M(E) = 0 .
   

   Now, the task

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

   Investigate the convergence of your solution towards the exact result with respect to the r_max and r_min parameters (separately) as well as with respect to the parameters acc and eps of your ODE integrator.

3. (1 point) Quadratic interpolation line-search

   Implement the quadratic interpolation line-search as described in the book.