AlephZero said:
Sorry, but I can't read the paper in your link without paying for it.
It's just a generalization of CN with effect of:
* having L-diagonal matrices (rather than 3-diagonal for CN)
* having s time-steps iterations (rather than 1 for CN)
The later form they put as repeated applications of the LU, so you get psi_{n+i/s} for i =1, ..., s, as intermediate steps. By computing the propagator once, you step directly to psi_{n+1}. I fail to see the advantage of the first case so I came here and asked.
AlephZero said:
If the answers are the same to 13 figures, the maths of the algorithm must be right (or you got VERY liucky!)
I don't think it's luck because it works with different potentials, different boundary conditions, different wavefunctions, etc.
AlephZero said:
Your description of the Mathematica solver seems to be complete overkill for this. You shouldn't need any pivoting, iterative improvement, etc.
That's the method it uses, it's based on LU, because that's what is most efficient. I copied this bit so that the answer doesn't come out to be "Mathematica just use completely inadequate algorithms in one case". It's using LU.
AlephZero said:
Try using the code from a library like LINPACK or LAPACK, and pick the routines that are the best match to the form of your equations (i.e. if you have banded symmetric positive definite equations, don't use the general purpose full matrix solver).
I'm not trying to get the thing working, I've got it running to the accuracy I want. The problem I have is that I see a complete contradiction with the literature that insists on doing it in a given way, by repeated applications of the LU decompositions, whereas a naive approach (easier to write down) turns out to be much more efficient. So I am puzzled why people insists on one way and not the other. You tell me it's not numerical accuracy (which I don't see as a problem either) and that in principle, from the complexity of the algorithm, LU is orders of magnitude better than inverting the matrix.
The only logical explanation is that Mathematica turns out to be very inefficient to do the LU (which I find very spooky, their linear algorithm are really good, they automatically detect sparse matrices and so on) and if I would do it in low-level language, I would go from 140s to 0.001s and get better time with LU as compared to inversion of matrix which would stay at 0.1s or so.
I'll have to take that on faith, or test more with low-level code like c or python. Still it's perplexing that LU is systematically put as *the* way to do it, when in some cases at least, it turns out to be orders of magnitude slower for the same accuracy. If I am correct, which I won't take as granted yet, people should say: "in some cases, like if you do this with a high-level language like Mathematica, it turns out to be better to compute the propagator once and merely iterate it rather than going through a series of LU decompositions".