# csape
带端点条件的三次样条插值
函数库: TyCurveFitting
# 语法
pp = csape(x,y)
pp = csape(x,y,conds)
pp = csape(x,[e1,y,e2],conds)
pp = csape((x1,...,xn),___)
# 说明
pp = csape(x,y) 以 pp 型返回给定数据 (x,y) 的三次样条插值。该函数将拉格朗日端点条件应用于数据的每一端,并将样条端点斜率与拟合每端最后四个数据点的三次多项式的斜率相匹配。位于同一位点的数据值将会取平均值。
pp = csape(x,y,conds)使用由 conds 指定的端点条件。示例
pp = csape(x,e1,y,e2,conds)使用由 conds 指定的端点条件,值e1和e2分别用于匹配左右端点值。示例
pp = csape((x1,...,xn),___) 使用单变量网格输入 x1,...,xn 返回网格数据的三次样条插值。在这种情况下,y 是一个 n+r 维数组,其中 r 是每个数据值的维度。conds 是一个带有 n 个元素的元组,它为每个 n 变量提供端点条件。在某些情况下,必须为端点条件提供端点条件。可以将此语法与先前语法中的任何参数一起使用。示例
# 示例
使用csape自定义端点条件
可以使用 csape 函数实现自定义端点条件。假设想在最左边的端点强制执行以下条件,e = x[1]
对于给定的标量 a、b 和 c。此时可以使用
在
此示例使用钛测试数据,这是数据拟合中使用的标准数据集。使用 titanium 函数加载数据。
using TyCurveFitting
using TyPlot
x,y = titanium()
定义 λ(s) 的系数。
a = -2
b = -1
c = 0
端点条件适用于数据集的最左端。
e = x[1]
现在,计算数据集的三次样条插值而不施加边界条件。
s1 = csape(x,y)
为计算
yZero = zeros(1,length(y))
1×2矩阵 conds 通过指定要固定的样条导数来设置边界条件。此示例仅需在数据的左端使用边界条件,因此这里用 conds 固定最左端的一阶导数。在右端,则固定函数本身的值。
conds = [1 0]
要指定将函数或其导数固定到的值,请将它们作为附加值添加到数据集参与拟合 - 在本例中为 yZero。 第一个元素指定左端的值,而最后一个元素指定右端的值。
在左端,将样条的一阶导数固定为 1。在右端,将函数值固定为 0(yZero 的最后一个元素的原始值)。在 yZero 的相应末端连接上这些边界条件值,并使用 csape 寻找满足这些边界条件且匹配数据值的样条。
s0 = csape(x,[1 yZero 0],conds)
使用之前提的 s 的表达式计算完全符合条件的样条。为了做到这点,使用样条
d1s1 = fnder(fnbrk(s1, 1))
d2s1 = fnder(d1s1)
d1s0 = fnder(fnbrk(s0, 1))
d2s0 = fnder(d1s0)
计算样条第一个多项式片段的导数,因为端点条件仅对于最左边的数据。
lam1 = a * fnval(d1s1, e)[1] + b * fnval(d2s1, e)[1]
lam0 = a * fnval(d1s0, e)[1] + b * fnval(d2s0, e)[1]
现在使用
pp = fncmb(s0, (c - lam1) / lam0, s1)
绘制样条,比较默认拟合与现有端点条件的结果。
fnplt(pp, [594, 632])
hold("on")
fnplt(s1, "b--", [594, 632])
plot(x, y, "ro", markerfacecolor = "r")
hold("off")
axis([594, 632, 0.62, 0.655])
legend(["Desired End Conditions", "Default end-conditions", "Data"], loc="southeast")
在第一个数据点附近的稳定点说明端点条件成功得到满足。
拟合多元数据
使用 csape 拟合多变量向量值数据。此示例对每个自变量使用不同的端点条件来拟合向量值数据。
首先,定义数据。对于此示例,在 2 维场中定义 3 维向量 v,在 x 方向具有固定条件或有规定的斜率,在 y 方向具有周期性端点条件。
using TyCurveFitting
using TyPlot
x = 0:4
y = -2:2
s2 = 1/sqrt(2)
v = zeros(3,7,5)
v[1,:,:] = [1 0 s2 1 s2 0 -1]'*[1 0 -1 0 1]
v[2,:,:] = [1 0 s2 1 s2 0 -1]'*[0 1 0 -1 0]
v[3,:,:] = [0 1 s2 0 -s2 -1 0]'*[1 1 1 1 1]
v 是一个3维数组,其中 v[:,i+1,j] 为坐标 x[i],y[j] 处的向量值。x 中的两个额外元素指定了斜率值:数据点 v[:,1,j] 和 v[:,7,j] 为固定端点条件提供沿线 x = 0 和 x = 4 的一阶导数的值。在 y 中,周期边界条件不需要任何额外的点。
现在,使用 csape 计算多元三次样条插值。
sph = csape((x,y),v,("clamped","periodic"))
要绘制结果,首先在合适的区间内计算样条。
val = fnval(sph,(0:.1:4,-2:.1:2))
surf(val[1,:,:],val[2,:,:],val[3,:,:])
daspect([1,1,1])
axis("off")
还可以使用 fnplt(sph) 计算和绘制样条曲面。注意,v 是一个 3 维数组,并且 v[:,i+1,j] 是在 (x[i],y[j]), i=1:5,j=1:5 处匹配的 3 向量。此外,因为 conds[1] 是"clamped",size(v,2) 是 7(而不是 5),v[r,:,j] 的第一和最后一个元素指定端点斜率值。
为端点条件提供端点条件
在某些情况下,需要为端点条件提供端点条件。在这个二变元的例子中,可以使用完全的双三次插值再现双三次多项式
using TyCurveFitting
using TyBase
sites = ([0, 1], [0, 2])
coefs = zeros(4, 4)
coefs[1, 1] = 1
g = ppmak(sites, coefs)
dorder = [1, 1]
Dxg = fnval(fnder(g, [1, 0]), sites)
Dyg = fnval(fnder(g, [0, 1]), sites)
Dxyg = fnval(fnder(g, [1, 1]), sites)
f = csape(sites, [Dxyg[1, 1] Dxg[1, :]' Dxyg[1, 2]; Dyg[:, 1] fnval(g, sites) Dyg[:, 2]; Dxyg[2, 1] Dxg[2, :]' Dxyg[2, 2]], ("complete", "complete"))
h = fnbrk(f, "c")
if any(squeeze(h) - coefs .!= 0)
println("something went wrong")
end
什么都没有打印出来。
# 输入参数
x — 数据位点向量 | 数组元组
要拟合的数据值 y 的数据位点,单变元时指定为一维数组,多变元时指定为元组。该函数在每个数据位点 x 创建样条 s 的节点,使得对于所有的 j,有 s[x[j]] = y[:,j]。
对于多元网格数据,指定 x 为元组可以为每个变量维度提供数据位点,使得 s[x1[i],x2[j],...,xn[k]] = y[:,i,j,...,k]。
数据类型: Int | Float
y-拟合数据向量 | 矩阵 | 多维数组
创建样条时拟合的数据值,指定为向量、矩阵或多维数组。可以将数据值 y[:,j]指定为标量、矩阵或 n 维数组。在同一数据位点 x 给出的数据值会被平均。
数据类型: Int | Float
conds - 端点条件"clamped" | "complete" | "not-a-knot" | "periodic" | "second" | "variational" | 1×2 矩阵
样条端点条件,指定为 "clamped","complete","not-a-knot","periodic","second","variational",或者1×2 矩阵。conds 的预定义选项为数据的每端都加上相同的端点条件。可以提供 1×2 矩阵为两端提供不同的端点条件。
可用的预定义端点条件如下
| "clamped" 或 "complete" | 将端点斜率与给定的值 e1 和 e2 匹配。如果不提供 e1 和 e2 的值,则此选项匹配默认的拉格朗日边界条件。 |
| "not-a-knot" | 忽视第二个和倒数第二个位点。此选项会忽略 e1 和 e2 被提供的任何值。 |
| "periodic" | 将左端的一阶和二阶导数与右端的一阶导数和二阶导数匹配。 |
| "second" | 将端点的二阶导数与给定的值 e1 和 e2 匹配。如果不提供 e1 和 e2 的值,则此选项对两者都使用默认值 0。这种情况下,此选项与"variational"相同。 |
| "variational" | 将端点二阶导数设置为零。此选项会忽略 e1 和 e2 被提供的任何值。 |
要在每一端指定不同的端点条件,将 conds 以 1×2 矩阵形式提供。此矩阵元素的元素指定由边界条件固定的样条导数的阶数。设置 conds[j] = i 将第 i 个导数
默认端点条件值是左四个位置的三次插值的导数,当 conds[1] = 1 时,其余情况则是 0。通过分别指定 e1 和 e2 来设置数据左端和右端的端点条件值。
可以指定 conds[j]为 0,1 或 2。如果指定不同的值或不指定 conds[j],则 conds[j] 是1,并且相应的边界条件值是默认值。
可以使用以下预定义的端点条件
| "clamped" | 如果 conds[j] == 1 | |
| "curved" | 如果 conds[j] == 1 | |
| "periodic" | 如果 conds = [0 0] | |
| "variational" | 如果 conds[j] == 2 && ej == 0 |
e、a 和 b 指的是数据左端或右端位置;ej 在数据的左端取 e1,在右端取 e2。
可以提供可选的端点条件值 e1 以及 e2,无论是否使用预定义或用户定义的选项 conds。但是,注意 conds 的一些预定义选项会忽略提供的端点边界条件值。
例如: "clamped" ,[1 0]
数据类型: String
e1 — 左端条件值标量
样条的左端条件值,指定为标量。e1 指定数据左端的第 i 个导数的值,其中conds提供 i。即使在两端使用不同的端点条件,如果在一端提供端点条件值,也必须为另一端提供一个端点条件值。
注意,一些 conds 的预定义选项会忽略所提供的任何边界条件值。
当 conds[1] = 1时,e1 的默认值是左四个位置的三次插值的导数,否则为 0。
数据类型: Int | Float
e2 — 右端条件值标量
样条的右端条件值,指定为标量。e2 指定数据右端的第 i 个导数的值,其中 conds 提供 i。即使在两端使用不同的端点条件,如果在一端提供端点条件值,也必须为另一端提供一个端点条件值。
请注意,一些 conds的预定义选项会忽略您提供的任何边界条件值。
当 conds[2] = 1时,e2 的默认值是右四个位置的三次插值的导数,否则为 0。
数据类型: Int | Float
# 输出参数
pp-样条结构体样条结构体
pp型的样条,为带有以下元素的结构体。关于 pp 样条的更多信息,参照 pp 型。
Form - 样条格式"pp"
样条的格式,以 "pp" 返回。"pp" 样条是以分段多项式的形式给出。
数据类型: String
Breaks-样条节点位置向量 | 数组元组
样条节点位置,单变元时以向量,多变元时以向量元组形式返回。向量包含严格增元素,逐点分别表示各段多项式定义区间的左右端点。
数据类型: Int | Float
Coefs - 多项式的系数矩阵 | 数组
各段多项式的系数,单变元时以矩阵,多变元时以数组形式返回。
数据类型: Int | Float
Pieces - 多项式分段段数标量 | 向量
描述样条的多项式分段的段数,单变元时以标量,多变元时以各变元段数组成的向量形式返回。
数据类型: Int
Order - 多项式阶数标量 | 向量
描述样条的分段多项式的阶数,单变元时以标量,多变元时以各变元阶数组成的向量形式返回。
数据类型: Int
Dim - 维数标量
目标函数的维数,以标量形式返回。
数据类型: Int
# 算法
使用 Syslab 的稀疏矩阵相关功能构造和求解对应的三对角矩阵。
csape 指令为 PGS 进程 CUBSPL 的大幅扩张版。