边界处有约束的样条

Question

我在三维网格上测量了数据，例如f(x, y, t)。我想用样条在 t 的方向上插值和平滑这些数据。目前，我使用 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 处具有消失的导数。有没有办法强制执行此约束？

Answer 1

你的例子不工作（在python 2.7.9），所以我只勾勒出我的想法：

计算sp
通过 sp.derivative 求导并在相关时间（可能与您测量数据的时间相同）对其进行评估
将相关点设置为零（例如t=0时的值）
根据导数值计算另一个样条。
集成样条函数。我想您将不得不以数字方式执行此操作，但这应该不是问题。不要忘记添加一个常量，以获得您的原始功能。

Answer 2

我能想到的最好的办法是使用 scipy.optimize.minimize 进行约束最小化。样条的导数很容易，所以约束很简单。我会使用常规样条拟合 (UnivariateSpline) 来获得结点 (t)，并保持结点固定（当然还有度数 k），并改变系数 c。也许还有一种方法可以改变结的位置，但我会把它留给你。

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

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

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

Answer 3

这是一种方法。基本思想是用 splrep 获取样条的系数，然后在调用 splev 之前修改它们。样条曲线中的前几个节点对应于 x 值范围内的最低值。如果对应于它们的系数设置为彼此相等，则在该端完全拉平样条曲线。

使用与示例中相同的数据、时间、x、y：

# set up example data
data = np.arange(1000).reshape((5, 5, 40))
times = np.arange(data.shape[-1])
y = 3
x = 3

# make 1D spline
import scipy.interpolate
from pylab import * # for plotting
knots, coefficients, degree = scipy.interpolate.splrep(times, data[y, x])
t = linspace(0,3,100)
plot( t, scipy.interpolate.splev(t, (knots, coefficients, degree)) )

# flatten out the beginning
coefficients[:2] = coefficients[0]
plot( t, scipy.interpolate.splev(t, (knots, coefficients, degree)) )
scatter( times, data[y, x] )
xlim(0,3)
ylim(720,723)

蓝色：原始点和通过它们的样条曲线。绿色：修改后的样条曲线，开头导数=0。两者都放大到最开始。

plot( t, scipy.interpolate.splev(t, (knots, coefficients, degree), der=1), 'g' )
xlim(0,3)

调用splev(..., der=1)绘制一阶导数。导数从零开始并稍微超调，以便修改后的样条可以赶上（这是不可避免的）。

修改后的样条曲线没有通过它所基于的前两个点（它仍然准确地命中所有其他点）。可以通过在原点旁边添加一个额外的内部控制点来修改它，以获得零导数并通过原始点；试验节点和系数，直到它达到你想要的效果。

边界处有约束的样条

Spline with constraints at border

python

scipy

spline