Re: interpolate

Nils Wagner wrote:
> Hi all,
> 
> Let us consider a curve which is given by m points (x_1,y_1) ...
> (x_m,y_m).
> How can I find the first and second derivative of this curve in
> x_1,...,x_m ?
> A small example illustrating the solution to this task would be
> appreciated.
> 
> Thanks in advance
> 
>                   Nils

I'd suggest using Neville's algorithm, which allows you to fit an exact 
polynomial of degree N-1 to a set of N points.  Finding the first and 
second derivatives should then be straight forward.  Usually it isn't 
necessary to fit all N points, but to take 2 (3) points on either side 
of the region that you want to interpolate and then fit the 4 (6) points 
using a 3rd (5th) order polynomial.  I've attached a Python version of 
the algorithm.

I learned about this algorithm a couple of years ago and have wondered 
how I did with out for so long.  This is real forgotten gem as far as 
I'm concerned.

Cheers,
Paul

-- 
Paul Barrett, PhD      Space Telescope Science Institute
Phone: 410-338-4475    ESS/Science Software Group
FAX:   410-338-4767    Baltimore, MD 21218

def neville(x, y, x0):
    """Through any N points y[i] = f(x[i]), there is a unique
    polynomial P order N-1.  Neville's algorithm is used for finding
    interpolates of this unique polynomial at any point x."""

    n = len(x)
    p = n*[0]
    for k in range(n):
        for j in range(n-k):
            if k == 0:
                p[j] = y[j]
	    else:
                p[j] = ((x0-x[j+k])*p[j]+(x[j]-x0)*p[j+1])/(x[j]-x[j+k])
    return p[0]

xa = [4,  12,  20]
print 'xa =', xa
xi = [0,   8,  16,  24]
print 'xi =', xi

l1, l2 = 8, 16
for m in range(0, 10, 1):
    #ya = [0.,  1., 1.*m]
    #print 'ya =', ya,
    #z1 = neville(xa, ya, l1)
    #z2 = neville(xa, ya, l2)
    #print l1, z1, l2, z2, (z2-z1), (ya[1]-(z2-z1)/2.)/2.,(ya[1]+(z2-z1)/2.)/2.

    yi = [0.,  0., 1., 1.*m+1.]
    print 'yi =', yi,
    z = neville(xi, yi, 12)
    s1, s2 = 2*(z-yi[1]), 2*(yi[2]-z)
    d = s2-s1
    print z, s1, s2, d, s1-d/2.

#m, d = 100, 4
#y = [0., 0.,  m*1.,  (m+1)*1.]
#y = [0., 24.9375,  25.9375,  25.9375]
#y = [0., 8.0430,  9.0430,  9.0430]
#y = [0., 2.7634,  3.7638,  3.7638]
#for l in range(0, 24, d):
#    print l+d, neville(x, y, l+d)-neville(x, y, l)