Numerical Recipies Spline

Question

I am trying to implement the numerical receipes function spline and splint as described here:

http://www.arcetri.astro.it/irlab/library/recipes/bookcpdf/c3-3.pdf -- Section 3.3

The code is still pretty confusing to me. Can someone explain what this number means?

if (yp1 > 0.99e30) -- What is this number?

Or if anyone could give some links to explain this more, I would be grateful.

Thanks.

Answer 1

Yp1 and Ypn is the first derivitive of the interpolating function at 1 and n respectivily. so derivitive of interpolating function at 1 = Yp1

From looking at the link again this is my opinion on why they use (yp > 0.99e30) The value of u[] is stored as a float. The smallest value that can be stored in a float is 1.2e-38(according to this link http://floating-point-gui.de/formats/fp/) . Now the function for calculating u[1] is

u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
Let (3.0/(x[2]-x[1]) = A, (y[2]-y[1]) = B, (x[2]-x[1])-yp1) = C.
u[1]= A * (B/C);

Now if yp is a very large number (eg yp > 1e30) This means C is a very large number. Which means B/C is a very small number which means A*(B/C) ~= 0(to small to be stored accurately); So there basically saying if the first derivative is too large(or too small) to be used where not going to use it (which is called natural spline and has an error of O(h^4)) and if they can use it they are going to use a clamped spline (which has an error of O(h^2) which means its more accurate than the natural spline method)

The same goes for Ypn.

Note you should add in another check to if(Yp > 10e30) to check for negative values. such as.

 if((Yp > 10e30 || Yp < -10e30))

Same for Ypn