[SciPy-User] Forced derivative interpolation??

Joe Kington jkington@wisc....
Fri Oct 9 13:40:47 CDT 2009


Hi Marco,

First off, take a close look at what Anne recommended.  It's probably much
more efficient to do this using another method if there is one.

Also, as Robert said, please reply using the title of your original post,
not the title of the digest.  It makes things much easier to find! :)

Keep in mind that there are many ways to do exactly what you've described.

As with any interpolation method, there are an infinite number of "right"
ways to do it that will all yield different results.

The simplest way to do this is to fit a piecewise cubic function to each
pair of points and derivatives.  (In other words, a cubic spline, but that
can mean many things)

So what we have is position (x) as a function of time(t):
x(t) = At^3 + Bt^2 + Ct + D
For each pair of points(t_0, x_0, t_1, x_1) and their derivatives (v_0,
v_1), we want to solve for A,B,C,&D such that:
x_0 = At_0^3 + Bt_0^2 + Ct_0 + D
x_1 = At_1^3 + Bt_1^2 + Ct_1 + D
v_0 = 3At_0^3 + 2Bt_0 + C
v_1 = 3At_1^3 + 2Bt_1 + C

Now, I'm lazy and don't particularly feel like doing the algebra, so let's
just put this in matrix form:
|x_0|      | t_0^3  t_0^2  t_0 0|   |A|
|x_1|  =   | t_1^3  t_1^2  t_1 0| * |B|
|v_0|      |3t_0^2  2t_0   1   0|   |C|
|v_1|      |3t_1^2  2t_1   1   0|   |D|
Which we can express as d = G*m and solve with numpy.

Now this isn't particularly elegant code, and I apologize for the lack of
comments (and general lack of clarity), but it should work as an example.

See the attatched image for an example plot.

import scipy as sp
import numpy as np
from matplotlib import pyplot as plt

def main():
    t = np.arange(11)
    x = np.random.randn(t.size)
    v = np.random.randn(t.size)
    ti = np.linspace(-1,11,100)
    xi = interpolate(t, x, v, ti)
    plot(ti, xi, t, x, v)

def interpolate(t, x, v, ti):
    """Given a set of points and the expected derivatives at those points,
    fit a cubic spline using the points and derivaties as constraints
    and return the interpolated values for each point in "ti"
    Input:
        t: The independent variable at each observed point
        x: The dependent variable at each observed point
        v: The derivate dx/dt at each observed point
        ti: The (new) t-values to solve for x at
    Output:
        xi: The interpolated x-values
    """
    # Note: there's probably a more elegant way of doing this...
    xi = np.zeros(ti.shape)
    for i in range(1, len(t)):
        t0, t1 = t[i-1], t[i]
        G = np.array([
            [t0**3, t0**2, t0, 1],
            [t1**3, t1**2, t1, 1],
            [3 * t0**2, 2 * t0, 1, 0],
            [3 * t1**2, 2 * t1, 1, 0],
            ])
        d = np.array([x[i-1], x[i], v[i-1], v[i]])
        m = np.linalg.solve(G,d)
        # One way of setting boundary conditions...
        if i == 1: # Points to the left of our data range
            interval = ti <= t[i]
        elif i == len(t)-1: # Points to the right of our data range
            interval = ti >= t[i-1]
        else: # Within bounds of data
            interval = (ti >= t[i-1]) & (ti <= t[i])
        t_int = ti[interval]
        xi[interval] = m[0] * t_int**3 + m[1] * t_int**2 + m[2] * t_int +
m[3]
    return xi

def plot(x,y, x0, y0, velocities):
    plt.figure()
    plt.hold(True)
    plt.plot(x, y, 'b-')
    plot_slopes(x0, y0, velocities)
    plt.plot(x0, y0, 'ro')
    plt.title('Fitting a spline using derivatives as constraints')
    plt.xlabel('Time')
    plt.ylabel('Position')
    plt.show()

def plot_slopes(x, y, velocities):
    from math import atan, cos, sin
    marker_width = np.diff(x).mean() / 4
    for x0, y0, v in zip(x,y,velocities):
        xdist = marker_width * cos(atan(v))
        ydist = marker_width * sin(atan(v))
        xleft, xright = x0 - xdist, x0 + xdist
        yleft, yright = y0 - ydist, y0 + ydist
        plt.plot([xleft, xright], [yleft, yright], 'r-')

main()


On Wed, Oct 7, 2009 at 2:44 AM, Marco Nicoletti <
nicoletti@consorzio-innova.it> wrote:

> Dear Joe,
>
> what I want to do is to interpolate the position array adding the
> constrain about the first derivative (the velocity array) whiche I
> already have.
>
> An example: I have t = [0,2,4,6] (the definition domain), p(t) =
> [10,13,16,18] (the position array) and v(t) = [1.1, 0.8, 0.7, 0.4].
> I have this new domain t1 = [0,1,2,3,4,5,6] and I want to obtain p1(t1)
> imposing the constrain about the velocities v1(t) = v(t).
> I want the spline interpolation.
>
> Any ideas?
>
> Thanks in advanced for your advices.
>
> Marco Nicoletti
>

On Wed, Sep 30, 2009 at 2:12 PM, Anne Archibald
<peridot.faceted@gmail.com>wrote:

> 2009/9/30 Marco Nicoletti <nicoletti@consorzio-innova.it>:
> > Dear all,
> >
> > I want to implement a spline interpolation forcing the condition on the
> > first or second derivative.
> > In other words I have a vector of position (p), velocity (v) and
> > acceleration (a) values;
> > I want to interpolate the position (p) vector imposing the conditions on
> the
> > velocity and acceleration values.
> >
> > The class UnivariateSpline() or intrp1D() in scipy.interpolate package
> don't
> > take as parameter the derivatives
> > (they export a method to evaluate derivatives).
> >
> > Any suggestions?
>
> If I have correctly understood your question, what you want to do is
> produce an interpolating spline with not just specified point values
> but specified derivative values at the given points. Scipy has at
> least two different pieces of code that might help. The first is, in
> recent versions of scipy, scipy.interpolate.PiecewisePolynomial. This
> allows you to fit a piecewise polynomial through a set of points,
> specifying derivatives at each point. It doesn't allow you to impose a
> spline-like constraint that higher derivatives must be continuous at
> the points. Its evaluation is also implemented in pure python, so it
> won't be terribly fast.
>
> A second option, useful if you need fast evaluation, is to abuse
> scipy's spline functions. scipy.interpolate.splrep doesn't take
> derivatives, but what it returns is a triple t, c, k. Given a t, c, k,
> you can then call splev, splint, splder, etcetera to get nice fast
> evaluation in compiled code. So what you can do is fabricate your own
> t, c, and k values. t is the list of knots, c is some sort of
> coefficients, and k is the order of the spline. The brute-force way I
> found to get these splines to produce the derivatives I wanted
> required me to repeat values in the t array. But once you've fixed the
> t array, the result is linear in the c values, so a little trial and
> error will give you formulas to produce any curve you need.
>
> Good luck,
> Anne
> > Thanks very much and have a nice day!
> >
> > Marco Nicoletti
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User@scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> >
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20091009/944cb562/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spline_derivative_example.png
Type: image/png
Size: 33686 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20091009/944cb562/attachment-0001.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: forced_derivative_interpolation.py
Type: text/x-python
Size: 2370 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20091009/944cb562/attachment-0001.py 


More information about the SciPy-User mailing list