边界有样条线 [英] Spline with constraints at border

问题描述

我已经在三维网格上测量了数据，例如f(x, y, t).我想用样条线在t方向上插值和平滑此数据. 目前，我使用scipy.interpolate.UnivariateSpline:

I have measured data on a three dimensional grid, e.g. f(x, y, t). I want to interpolate and smooth this data in the direction of t with splines. Currently, I do this with scipy.interpolate.UnivariateSpline:

import numpy as np
from scipy.interpolate import UnivariateSpline

# data is my measured data
# data.shape is (len(y), len(x), len(t))
data = np.arange(1000).reshape((5, 5, 40))  # just for demonstration
times = np.arange(data.shape[-1])
y = 3
x = 3
sp = UnivariateSpline(times, data[y, x], k=3, s=6)

但是，我需要花键在t=0处具有消失的导数.有没有办法强制执行此约束?

However, I need the spline to have vanishing derivatives at t=0. Is there a way to enforce this constraint?

推荐答案

我能想到的最好的事情是使用scipy.optimize.minimize在约束条件下进行最小化.取样条曲线的导数很容易，因此约束很简单.我将使用常规的样条拟合(UnivariateSpline)来获得结(t)，并保持结固定(当然，度数为k)，并更改系数c.也许还有一种方法可以改变结的位置，但我会把它留给您.

The best thing I can think of is to do a minimization with a constraint with scipy.optimize.minimize. It is pretty easy to take the derivative of a spline, so the constraint is simply. I would use a regular spline fit (UnivariateSpline) to get the knots (t), and hold the knots fixed (and degree k, of course), and vary the coefficients c. Maybe there is a way to vary the knot locations as well but I will leave that to you.

import numpy as np
from scipy.interpolate import UnivariateSpline, splev, splrep
from scipy.optimize import minimize

def guess(x, y, k, s, w=None):
    """Do an ordinary spline fit to provide knots"""
    return splrep(x, y, w, k=k, s=s)

def err(c, x, y, t, k, w=None):
    """The error function to minimize"""
    diff = y - splev(x, (t, c, k))
    if w is None:
        diff = np.einsum('...i,...i', diff, diff)
    else:
        diff = np.dot(diff*diff, w)
    return np.abs(diff)

def spline_neumann(x, y, k=3, s=0, w=None):
    t, c0, k = guess(x, y, k, s, w=w)
    x0 = x[0] # point at which zero slope is required
    con = {'type': 'eq',
           'fun': lambda c: splev(x0, (t, c, k), der=1),
           #'jac': lambda c: splev(x0, (t, c, k), der=2) # doesn't help, dunno why
           }
    opt = minimize(err, c0, (x, y, t, k, w), constraints=con)
    copt = opt.x
    return UnivariateSpline._from_tck((t, copt, k))

然后我们生成一些初始斜率为零的假数据并对其进行测试:

And then we generate some fake data that should have zero initial slope and test it:

import matplotlib.pyplot as plt

n = 10
x = np.linspace(0, 2*np.pi, n)
y0 = np.cos(x) # zero initial slope
std = 0.5
noise = np.random.normal(0, std, len(x))
y = y0 + noise
k = 3

sp0 = UnivariateSpline(x, y, k=k, s=n*std)
sp = spline_neumann(x, y, k, s=n*std)

plt.figure()
X = np.linspace(x.min(), x.max(), len(x)*10)
plt.plot(X, sp0(X), '-r', lw=1, label='guess')
plt.plot(X, sp(X), '-r', lw=2, label='spline')
plt.plot(X, sp.derivative()(X), '-g', label='slope')
plt.plot(x, y, 'ok', label='data')
plt.legend(loc='best')
plt.show()

这篇关于边界有样条线的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！