[SciPy-user] interpolate

Nils Wagner nwagner at mecha.uni-stuttgart.de
Fri Aug 1 09:35:58 EDT 2003


Paul Barrett schrieb:
> 
> 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
> 
>   ------------------------------------------------------------------------
>                  Name: neville.py
>    neville.py    Type: Plain Text (text/plain)
>              Encoding: 7bit

Thank you for your reply. Finally, I have used interpolate.splrep and
interpolate.splev.
I have enclosed a short demo:

from scipy import *
from scipy.xplt import *
import gui_thread
x = arange(0,1,0.1)
#
# Original function
#
y = exp(-x/3.0)
#
# Analytical derivatives
#
ys  = -exp(-x/3.0)/3.0
yss = +exp(-x/3.0)/9.0
#
tck = interpolate.splrep(x,y)
#
# Spline
#
z=interpolate.splev(x,tck,0)
#
# First derivative
#
zs=interpolate.splev(x,tck,1)
#
# Second derivative
#
zss=interpolate.splev(x,tck,2)

window(0)
xplt.plot(x,y,'x',x,z,'ro')
window(1)
xplt.plot(x,ys,'x',x,zs,'ro')
window(2)
xplt.plot(x,yss,'x',x,zss,'ro')

Nils




More information about the SciPy-User mailing list