8.  Boundary Value Problems -- matrix methods for linear systems

An alternative to shooting is simply to write the differential equation as a difference equation on a grid, and solve the entire system, boundary conditions included, as a matrix problem. For example, if our equation is

	    y'' + g(x) y = 0
	    y(a) = ya, y(b) = yb

we can difference on a uniform grid in x with N+1 points (indices 0 to N, N intervals) and spacing h = (b-a)/N:

	    h2yn'' = yn+1 - 2 yn + yn-1

and recast the system as a matrix equation of dimension N-1 involving the interior points (n = 1, ..., N-1):

	    A y = r

where the derivative matrix is

	    A = (-2+k1   1     0     0     0     0 ...)
	        (  1   -2+k2   1     0     0     0 ...)
[	        (  0     1   -2+k3   1     0     0 ...)
   	        (  0     0     1   -2+k4   1     0 ...)
	        (            ... etc. ....            )

where k_n = h² g(x_n) and

	    r = (-ya)
	        ( 0 )
	        ( 0 )
	        ( . )
	        ( . )
	        ( . )
	        (-yb)

takes care of the boundary conditions. You can find a more general derivation of the procedure here.

This equation is easily solved, for example using the solve function in scipy.linalg:

	    y = solve(A, r)

Here is a schematic example indicating how the pieces fit together.

• Exercise 8.1: Use the matrix method just described to solve the simple problem addressed in Exercise 7.1:

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

  Use a grid with 101 points (100 intervals). Plot the numerical and analytic solutions, and determine the maximum error on the grid. Repeat using a grid of 1001 points. How does the error scale with grid spacing?   Solution.

  Note that this method only applies to linear systems. We will see later how to approach nonlinear problems. Note also that since all the matrices we will encounter are tridiagonal, functions exist to solve them more efficiently, but for now we won't worry too much about speed.

• Exercise 8.2: Modify your program to include a centered approximation to the first derivative y', as described in more detail here, and solve the following system:

  	    y'' + y'/x + (1-1/x**2)y = 0
  	    y(0) = 0, y(10) = 1

  Use a grid with 1001 points (1000 intervals). In this case the solution should just be a multiple of the Bessel function J₁. Also solve the problem by shooting just to check!   Matrix solution.   Shooting solution.
  

  
  

• Application to eigenvalues

  The matrix method generalizes naturally to eigenvalue problems, so long as the boundary conditions are y(x) = 0 for x = a,b. This is somewhat restrictive in general, but as we have seen, it is almost always OK for the Schrödinger problem. For a system like

  	    y'' + g(x) y  + z y = 0
  	    y(a) = 0, y(b) = 0

  the problem reduces to finding the eivenvalues and eigenvectors of the matrix A, as defined above. The eigenvalues and eigenvectors can be determined using the function (from scipy.linalg) eig:

  	    eigenvalues,eigenvectors = eig(A)

  where eigenvectors[:,i] is the eigenvector corresponding to eigenvalue eigenvalues[i] (and, since the matrix is real and symmetric, there are exactly N-2 of each). Note that the eigenvalues are not sorted, so it is generally desirable to first sort using the numpy argsort function, which returns the indices of the eigenvalues (and eigenfunctions) sorted as desired. Although matrix diagonalization is an O(N³) process, the facts that the routines are pure numpy (i.e. C) and return all solutions at once can often make this approach competitive with the shooting algorithms discussed earlier.

• Exercise 8.3: Recast the matrix method described in the above examples to solve the eigenvalue problem

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

  on a grid with 101 points (100 intervals). Print out the lowest 4 eigenvalues you obtain and plot the corresponding eigenfunctions.   Here's a starting point. Solution.

• Exercise 8.4: Repeat Exercise 7.5, using matrix methods: 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

  
  

  
  

• The Numerov algorithm

  The matrix method just described is based on a second-order differencing of the second derivative:

  	    yn+1 - 2 yn + yn-1 = h2yn'' + O(h4)

  In 1924 a Russian astronomer named Numerov realised that, since the next term is in fact h⁴y''''/12, it can also be estimated by differentiating the original differential equation. Details can be found here. The result is a remarkably compact fourth-order method (albeit one only applicable to systems without a y' term, like the Schrödinger equation), that can be applied with essentially no additional work compared to the second-order scheme, but leads to substantially more accurate results.

  Here is an almost complete Numerov solver to use as a starting point.

• Exercise 8.5: Repeat Exercise 8.1 using the Numerov method.

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

  Use a grid with 11 points (10 intervals). Plot the numerical and analytic solutions, and determine the maximum error on the grid. Repeat using a grid of 101 points. How does the error scale with grid spacing?   Solution.

  With the Numerov method, the eigenvalue problem described above,

  	    y'' + g(x) y  + z y = 0
  	    y(a) = 0, y(b) = 0

  becomes the matrix equation

  	    A y + (h2z/12) B y = 0

  where

  	    A = ( β1  α1  0   0   0   0 ...)
  	        ( γ2  β2  α2  0   0   0 ...)
  	        ( 0   γ3  β3  α3  0   0 ...)
     	        ( 0   0   γ4  β4  α4  0 ...)
  	        (      ... etc. ....       )

  with  α_n = 1 + h²g_n+1/12, and  β_n = -2(1 - 5h²g_n/12),  γ_n = 1 + h²g_n-1/12, and

  	    B = ( 10   1   0   0   0   0 ...)
  	        (  1  10   1   0   0   0 ...)
  	        (  0   1  10   1   0   0 ...)
  	        (  0   0   1  10   1   0 ...)
  	        (       ... etc. ....       )

  See here for details. The problem thus reduces to finding the eigenvalues and eigenvectors of the matrix B^-1A.

• Exercise 8.6: Repeat Exercise 8.3 using the Numerov method. How many grid points are needed to ensure accuracy to at least three decimal places?

• Exercise 8.7: Repeat Exercise 8.4 using the Numerov method.
  
  
  
  

• Further development: brief discussion of

  • Basis expansion

  • Nonlinear systems

  • Context: Stellar structure