# numerov.R
#
# Use the Numerov algorithm to integrate the 2nd order differential
# equation,
#
#     P''(r) = Q(r)*P(r)
#
# The Numerov algorithm may be expressed by the recurrence
# relation,
#
#  F_n+1 = [(2 + 10*T_n)/(1 - T_n)]*F_n - F_n-1
#
# where,
#
#    T_n = (h^2/12)*Q_n
#
#    F_n = (1 - T_n)*P_n
#
# Given the input array, Q, and the first two values of the
# array P, P[1] and P[2], initialized, the recurrence relation
# is applied to compute the successive values of P[n]. 
#
#
# References:
#
# 1. B. Numerov, Publs. observatoire central astrophys. Russ., 
#    v. 2, p. 188 (1933).
#
# 2. http://en.wikipedia.org/wiki/Numerov%27s_method
#
# 3. Anders Blom, 2002, "Computing algorithms for solving the 
#    Schroedinger and Poisson equations", available on the web at
#    http://www.teorfys.lu.se/personal/Anders.Blom/useful/scr.pdf
#
# 4. E. Merzbacher, Quantum Mechanics, 3rd ed., Wiley (1998). 
#
# Copyright (c) 2010--2012 Krishna Myneni
#
# This code may be used for any purpose, as long as the
# copyright notice above is preserved.
#
# Revisions:
#   2012-10-08  km; ported from the Forth module, numerov.4th,
#                   and program numerov-test.4th.
#

numerov_integrate <- function( P, Q, h ) {
	h2_12 = h^2/12
	T_nm1 = h2_12*Q[1]
	T_n = h2_12*Q[2]
	F_nm1 = (1 - T_nm1)*P[1]
	F_n = (1 - T_n)*P[2]

	for (i in 3:length(Q)) {
	  T_np1 =  h2_12*Q[i]
	  F_np1 =  (2 + 10*T_n)*F_n/(1 - T_n) - F_nm1
	  P[i]  <<-  F_np1/(1 - T_np1)
	  T_nm1 =  T_n
	  T_n   =  T_np1
	  F_nm1 =  F_n
	  F_n   =  F_np1
	}	
}


# Test Numerov integration for the case of the H atom in the 1s state:
#
#       Q(r) = -2/r + 1
#
# Since, at small r, P(r) goes as r^(l+1) [4], for this example, the
# initial value of P and its first derivative may be initialized by
# setting P[1] and P[2] to the following values:
#
#  P[1] = r[1],  P[2] = r[1] + h
#
# This should give us the correct solution apart from a
# normalization constant. Proper normalization of the result is
# illustrated below.

## Set up arrays Q and P
h = 0.001
r = seq(1e-4, 10, by=h)
Q = -2/r + 1
P = array(0, c(length(Q)))
P[1] = r[1]
P[2] = P[1] + h

numerov_integrate(P, Q, h)

## Plot the result of the numerical integration, along with
## the actual radial function for the H atom 1s state.
P_1s <- function(r) { return( 2*r*exp(-r) ) }

plot(r, P_1s(r), type="l", col="green", lwd=4, xlab="r (a.u.)", ylab="P(r)")
lines(r, P, col="blue")

## Normalize the numerical integration result and overlay on plot
norm = sum(P*P)*h
P_norm = P/sqrt(norm)
lines(r, P_norm, lty=2, col="black")

title(main="Numerov Integration of the Radial Equation for H atom 1s State")
legend(x=6, y=0.7, legend=c("Analytic Sol.", "Numerov Int.", "Normalized Numerov"),
       col=c("green", "blue", "black"), lty=c(1,1,2), lwd=c(4,1,1))