jcline
12-14-2010, 05:51 PM
I don't understand how to assign indices for the matrix S(i,j) that describes dependence of the boundary conditions on the dependent variables in the relaxation routine. According to the description of the algorithm in the NR text, at the first boundary, indices are assigned as
S(j,n) = dE(j,1)/dy(n,1)
where n runs from 1 to N. (I am describing the fortran version; in the C version replace 1 by 0 here and M+1 by M below, but this doesn't have any bearing on the inconsistency to be described.) According to the conventions painstakingly explained in the text, this should actually read
S(j,n+N) = dE(j,1)/dy(n,1)
because derivatives with respect to y defined at the same lattice point as E go in the second half of the rectangular S(j,n) block. Now in the example routine that follows, this is what is actually done, so I was happy to think the omission of +N in the index was just a typo (though if so I would have thought this should be corrected by the 3rd edition; it's still there).
Now when we come to the second boundary condition, the text tells us that
S(j,n) = dE(j,M+1)/dy(n,M)
where n = 1 to N (as opposed to N+1 to 2N). So far so good; this agrees with the convention that derivatives with respect to the lattice point on the left, compared to the position index of E, go into the first (left) half of S(j,n). But in the example program, at the second boundary condition, n goes from N+1 to 2N, in contradiction to the description of the algorithm. Apparently the program works. What is the correct description, and how has anyone ever been able to use this horribly documented routine?
S(j,n) = dE(j,1)/dy(n,1)
where n runs from 1 to N. (I am describing the fortran version; in the C version replace 1 by 0 here and M+1 by M below, but this doesn't have any bearing on the inconsistency to be described.) According to the conventions painstakingly explained in the text, this should actually read
S(j,n+N) = dE(j,1)/dy(n,1)
because derivatives with respect to y defined at the same lattice point as E go in the second half of the rectangular S(j,n) block. Now in the example routine that follows, this is what is actually done, so I was happy to think the omission of +N in the index was just a typo (though if so I would have thought this should be corrected by the 3rd edition; it's still there).
Now when we come to the second boundary condition, the text tells us that
S(j,n) = dE(j,M+1)/dy(n,M)
where n = 1 to N (as opposed to N+1 to 2N). So far so good; this agrees with the convention that derivatives with respect to the lattice point on the left, compared to the position index of E, go into the first (left) half of S(j,n). But in the example program, at the second boundary condition, n goes from N+1 to 2N, in contradiction to the description of the algorithm. Apparently the program works. What is the correct description, and how has anyone ever been able to use this horribly documented routine?