Avoiding overflow in ``schint''

As we discussed in lab, we can exploit the linearity of the Schrödinger equation to avoid overflow of Psi and Psip in schint() by renormalizing the solution by multiplying it through by a small constant each time that the magnitude of the solution exceeds a critical value. The IEEE range of double precision values is somewhat beyond $10^{-300}$ to $10^{300}$. To ensure enough room to avoid overflows, it is safest to keep the magnitude of the solution between $10^{-150}$ and $10^{150}$. You may do this by scaling the solution by $10^{-300}$ every time its absolute value exceeds $10^{150}$. Once again, it is best to collect such special numerical constants in '#define' statements near to top of your file. You should thus add the following lines to the top of ``physics.c'' next to the '#define TOL (1.e-12)' statement,

#define SMALL (1.e-300)
#define SQRTSMALL (1.e-150)
#define BIG (1.e300)
#define SQRTBIG (1.e150)
Then, you should modify your main loop in schint(), so that whenever the absolute value of Psi (probably called y[1]) exceeds ``SQRTBIG'', your entire solution is multiplied by ``SMALL''.

Hint: When you rescale your solution, you must scale both components of $\vec y$ in your differential equation (i.e., both y[1] and y[2]), and you must also rescale all of the values previously generated and stored in Psiout[]. The rescaling of Psiout[] will require writing another loop and some care to be certain that all of the appropriate elements of Psiout[] are rescaled, for both dk$>$0 and dk$<$0.

Tomas Arias 2007-02-27