Python 中具有连续四阶和五阶导数的平滑样条？
Smoothing spline with continuous 4th and 5th order derivatives in Python?
我有一些嘈杂的时间序列，将输入参数表示为一组 ODE 的时间函数。这些嘈杂的时间序列需要进行插值，并且该插值函数的导数需要连续到与我的 ODE 求解器相同的阶数。否则，N 阶 ODE 求解器将尝试采用极小的自适应时间步来处理 N 阶导数中的不连续跳跃。因此，例如，对于 scipy.integrate.solve_ivp 中的标准 RK45 求解器，插值函数的所有导数至少达到 5 阶必须是连续的，我认为。
在 Python 中创建平滑样条以使其前 5 个导数连续的最佳方法是什么？我很难完成这个简单的任务。下面是一个最小的工作示例，显示即使使用来自 scipy 的 5 阶平滑 UnivariateSpline，4 阶导数显示 discontinuities/oscillations 并且 5 阶导数只是爆炸（~1e6 大小）。尽管我首先对嘈杂的时间序列进行高斯平滑以帮助插值函数，但还是如此。
我怀疑这至少是 scipy RK45 求解器花费如此长的时间来求解我的 ODE 系统的原因之一——它采用了不必要的大量较小的时间步长（尽管我认为我使用默认 solve_ivp 公差的事实也可能发挥作用）。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import lognorm
from scipy.ndimage import gaussian_filter1d
from scipy.interpolate import UnivariateSpline
x = np.sort(np.random.uniform(0.1,10,500))
dist = lognorm(0.5,loc=1.0)
y_noisy = dist.pdf(x) + np.random.uniform(0.0,1,500)
y_smooth = gaussian_filter1d(y_noisy,10,mode='nearest')
spl = UnivariateSpline(x,y_smooth,k=5,s=2)
xx = np.linspace(0.1,10,100000)
yy = spl(xx)
d1 = np.diff(yy) / np.diff(xx)
d2 = np.diff(d1) / np.diff(xx[1:])
d3 = np.diff(d2) / np.diff(xx[1:-1])
d4 = np.diff(d3) / np.diff(xx[1:-2])
d5 = np.diff(d4) / np.diff(xx[1:-3])
fig, axes = plt.subplots(nrows=6,ncols=1,figsize=(8,8))
fig.subplots_adjust(hspace=0.7)
axes[0].plot(x,y_noisy,'k-',lw=2,alpha=0.1)
axes[0].plot(x,y_smooth,'y-',lw=2)
axes[0].plot(xx, yy)
axes[0].set_title('5th-order smoothing UnivariateSpline')
axes[1].plot(xx[1:], d1)
axes[1].set_title('first derivative')
axes[2].plot(xx[1:-1], d2)
axes[2].set_title('second derivative')
axes[3].plot(xx[2:-1], d3)
axes[3].set_title('third derivative')
axes[4].plot(xx[3:-1], d4)
axes[4].set_title('fourth derivative')
axes[5].plot(xx[4:-1], d5)
axes[5].set_title('fifth derivative')
我认为这不是检查函数在五阶导数处是否连续的有效方法。我认为时间步太小了，四次和五次导数的图被浮点错误淹没了。第四个和第五个差异是如此之小，以至于小错误会产生大影响。
为了测试这一点，我用 np.sin() 替换了你拟合的样条，我知道它具有连续导数。
果然，它说四阶导数跳动了很多。但那不可能——正弦的四阶导数就是正弦。
如果我将时间步数减少到 1000，它会正确找到四阶和五阶导数。将相同的技术应用于原始函数，我得到：
我有一些嘈杂的时间序列，将输入参数表示为一组 ODE 的时间函数。这些嘈杂的时间序列需要进行插值，并且该插值函数的导数需要连续到与我的 ODE 求解器相同的阶数。否则，N 阶 ODE 求解器将尝试采用极小的自适应时间步来处理 N 阶导数中的不连续跳跃。因此，例如，对于 scipy.integrate.solve_ivp 中的标准 RK45 求解器，插值函数的所有导数至少达到 5 阶必须是连续的，我认为。
在 Python 中创建平滑样条以使其前 5 个导数连续的最佳方法是什么？我很难完成这个简单的任务。下面是一个最小的工作示例，显示即使使用来自 scipy 的 5 阶平滑 UnivariateSpline，4 阶导数显示 discontinuities/oscillations 并且 5 阶导数只是爆炸（~1e6 大小）。尽管我首先对嘈杂的时间序列进行高斯平滑以帮助插值函数，但还是如此。
我怀疑这至少是 scipy RK45 求解器花费如此长的时间来求解我的 ODE 系统的原因之一——它采用了不必要的大量较小的时间步长（尽管我认为我使用默认 solve_ivp 公差的事实也可能发挥作用）。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import lognorm
from scipy.ndimage import gaussian_filter1d
from scipy.interpolate import UnivariateSpline
x = np.sort(np.random.uniform(0.1,10,500))
dist = lognorm(0.5,loc=1.0)
y_noisy = dist.pdf(x) + np.random.uniform(0.0,1,500)
y_smooth = gaussian_filter1d(y_noisy,10,mode='nearest')
spl = UnivariateSpline(x,y_smooth,k=5,s=2)
xx = np.linspace(0.1,10,100000)
yy = spl(xx)
d1 = np.diff(yy) / np.diff(xx)
d2 = np.diff(d1) / np.diff(xx[1:])
d3 = np.diff(d2) / np.diff(xx[1:-1])
d4 = np.diff(d3) / np.diff(xx[1:-2])
d5 = np.diff(d4) / np.diff(xx[1:-3])
fig, axes = plt.subplots(nrows=6,ncols=1,figsize=(8,8))
fig.subplots_adjust(hspace=0.7)
axes[0].plot(x,y_noisy,'k-',lw=2,alpha=0.1)
axes[0].plot(x,y_smooth,'y-',lw=2)
axes[0].plot(xx, yy)
axes[0].set_title('5th-order smoothing UnivariateSpline')
axes[1].plot(xx[1:], d1)
axes[1].set_title('first derivative')
axes[2].plot(xx[1:-1], d2)
axes[2].set_title('second derivative')
axes[3].plot(xx[2:-1], d3)
axes[3].set_title('third derivative')
axes[4].plot(xx[3:-1], d4)
axes[4].set_title('fourth derivative')
axes[5].plot(xx[4:-1], d5)
axes[5].set_title('fifth derivative')
我认为这不是检查函数在五阶导数处是否连续的有效方法。我认为时间步太小了，四次和五次导数的图被浮点错误淹没了。第四个和第五个差异是如此之小，以至于小错误会产生大影响。
为了测试这一点，我用 np.sin() 替换了你拟合的样条，我知道它具有连续导数。
果然，它说四阶导数跳动了很多。但那不可能——正弦的四阶导数就是正弦。
如果我将时间步数减少到 1000，它会正确找到四阶和五阶导数。将相同的技术应用于原始函数，我得到：