7.  Boundary Value Problems -- shooting method

• Read Numerical Recipes sections 18.0 and 18.1.

• Introduction -- the Shooting method

• sample solver
• sample integrator

• Exercise 7.1: Combine these solvers to create a simple boundary-value problem solver. Apply it to the problem

  	    y'' = -y
  	    y(0) = 1, y(1) = 2

  What value of y'(0) solves the system? Add one or two applications of the secant method once bisection has converged to the desired tolerance to improve the accuracy of the method. Solution.  Animated solution.

• Exercise 7.2: Find a solution of

  	    y'' = -y - 10 y3
  	    y(0) = 1, y(1) = 2

  with 0 < y'(0) < 50. Be careful with your starting point for the bisection search. Solution

• Exercise 7.3: Find all solutions of

  	    y'' = -y - 10 y3
  	    y(0) = 1, y(1) = 2

  with -100 < y'(0) < 100, using the shooting method. Plot the solutions on a single graph, with a legend (or labels) stating which is which.

• For a linear eigenvalue equation, we typically have some boundary constraint on a combination of y and y', but no other global constraint on the solution -- that is taken care of by the overall normalization. Instead, the eigenvalue itself becomes the free parameter z we vary to satisfy the second boundary condition. The cleanest way to do this within the framework we have just developed is to add the eigenvalue z as an extra variable. That is, in the two dimensional case, we set y[2] = z and the system of equations to solve becomes (say)

  	    dy[0]/dx = y[1]
  	    dy[1]/dx = -y[2]*y[0]
  	    dy[2]/dx = 0

  Thus the function f must be modified to return a zero in its new (third) component. We continue to bisect on z using the (unchanged) function g(z). The function integrate changes from

              def integrate(func, a, ya, za, b, dx):
  	        x = a
  	        y = np.array([ya, za])
  	        while x < b-0.5*dx:
  	            x, y = rk4_step(f, x, y, dx)
  		return x, y

  to

              def integrate(func, a, ya, z, b, dx):
                  x = a
                  y = np.array([ya, 1.0, z])		# <---
                  while x < b-0.5*dx:
  		    x, y = rk4_step(f, x, y, dx)
                  return x, y

• Exercise 7.4: Find the smallest positive eigenvalue z for the equation

  	    y'' = -z y
  	    y(-1) = 0, y(1) = 0

  Solution  Animated solution.

• Exercise 7.5: Find the first 4 eigenvalues E and eigenfunctions y of the "confined harmonic oscillator" with

  	    y'' = -E y + U(x) y
  	    U(x) = x2, -5 < x < 5
  	    y(-5) = 0, y(5) = 0

  Solution

• What if the domain extends to infinity? In many cases, we can write down the solution at large x (or an approximation to it), and this can then determine the boundary condition to use. For example, in the finite square well problem, with

  	    U(x) = 0,  |x| < a
  	    U(x) = U0, |x| > a,

  and 0 < E < U0, we know that the exterior solution must be of the form

  	    y = A exp(-eta*x)

  where eta^2 = U0 - E = U0 - z here. Differentiating, we have

  	    y' + eta y = 0

  This is the boundary condition to be applied at |x| = a. This general approach applies to any system in which U(x) tends to a constant as |x| goes to infinity.

• A more complete treatment of the quantum finite well problem

• Exercise 7.6: Find an eigenvalue E and eigenfunction y of the finite square well with a = 1 and U0 = 4.

• Exercise 7.7: Find the ground state energy of a quantum system with
  applying the above boundary condition at |x| = 5.

• If |U(x)| is increasing as x goes to infinity, there is no simple boundary relation between y and y' and it is usually acceptable to set y = 0 at some large |x|. For example, in the case of the harmonic oscillator, it is easy to show that, for large |x|

              y  ~  e-x2/2

  and

              y'  ~  xy

  so we can take y = 0 at x = x0 (i.e. we simply adopt the "confined" oscillator presented earlier). In this case, however, the rapidly falling functional form of the solution means that shooting from -x0 to x0 can become numerically unstable, and it is better to shoot from the center to x = x0. This is particularly simple in cases (like all those we have seen) where the potential is symmetric about x = 0. We then expect even or odd solutions for y, which we can select by choosing

  	    y = 1,  y' = 0  (even)
  
  	    y = 0,  y' = 1  (odd)

  at x = 0.

• Summary of shooting strategies

• Exercise 7.8: Find the first 4 eigenvalues E and eigenfunctions y (2 even, 2 odd) of the harmonic oscillator with

  	    U(x) = x2

  by shooting from x = 0 to x = x0 = 6, with y(x0) = 0.

• Exercise 7.9: Find the even eigenfunction of the finite well problem (Exercise 7.6) by shooting from x = 0 and assuming y(x) = 0 at x = 10. Compare the eigenvalue you obtain with the result of Exercise 7.6.