1 Aug 2003 15:17

## 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.
>
>
>                   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.

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)
```

Gmane