利用有理样条进行保形插值是几何造型领域中的研究热点之一^[1-8]。文献[4-6]介绍的保形有理四次样条插值函数节点处导数的算法为构造光顺的保形有理四次样条插值曲线，以形状控制参数和节点处的导数为决策变量，以插值曲线应变能最小为目标函数，以插值函数保形以及形状控制参数和节点处的导数大于0作为约束条件，建立优化模型，给出插值算法是求解获得应变能最小的保形有理四次样条插值曲线。

1 插值函数^[4]

给定数据{(x_i, f_i), i=1, 2, …, n}，其中：f_i为被插值函数在节点x_i处的函数值，此处令x₁<x₂ < … < x_n。记h_i=x_i+1-x_i，Δ_i=(f_i+1-f_i)/h_i；令t=(x-x_i)/h_i，当x∈[x_i, x_i+1]时，t∈[0,1]，i=1, 2, …, n-1，构造分子是四次多项式、分母是一次多项式的有理样条插值曲线s(x)，s(x) 定义如下：

$ s\left( x \right) = s\left( {{x_i} + {h_i}t} \right) \equiv S\left( t \right)\left| {_{[{t_i}, {t_{i + 1}}]}} \right. = \frac{{{P_i}\left( t \right)}}{{{Q_i}\left( t \right)}}, i = 0, 1, 2, \ldots, n -1。 $

(1)

其中:

$ {P_i}\left( t \right) = {\alpha _i}{f_i}{\left( {1-t} \right)^4} + {U_i}{\left( {1-t} \right)^3}t + {V_i}{\left( {1-t} \right)^2}{t^2} + {W_i}\left( {1 - t} \right){t^3} + {\beta _i}{f_{i + 1}}{t^4}, $

(2)

$ \begin{array}{l} {Q_i}\left( t \right) = {\alpha _i}\left( {1-t} \right) + {\beta _i}t, {\rm{ }}\\ {U_i} = \left( {3{\alpha _i} + {\beta _i}} \right){f_i} + {\alpha _i}{d_i}{h_i}, {\rm{ }}\\ {V_i} = 3{\beta _i}{f_{i + 1}} + 3{\beta _i}{f_i}。 \end{array} $

(3)

插值函数s(x) 在节点x_i处的一阶导数值为d_i，α_i, β_i被称为形状控制参数且α_i>0, β_i>0。

由 (2)(3) 式易知，有理样条插值s(x) 满足下列插值性质：

$ s\left( {{x_i}} \right) = {f_i}, i = 0, 1, \cdots, n。 $

(4)

2 保单调性分析^[4-8]

如果f(x) 在区间[a, b]上有x₁ < x₂ < … < x_n, f₁≤f₂≤...≤f_n或Δ_i≥0，那么函数f(x) 在区间[a, b]上单调递增。选取导数值d_i时必须满足：

$ {d_i} \ge 0, i = 1, 2, \cdots, n。 $

(5)

在区间[x_i, x_i+1]上，s(x) 单调递增的充要条件为：

$ s\prime \left( x \right) \ge 0, i = 0, 1, \cdots, n-1。 $

(6)

对 (1) 式中定义的插值函数s(x) 求导：

$ s\prime \left( x \right) = S\prime \left( t \right){h_i}^{- 1} = \frac{{\sum\limits_{j = 0}^4 {{Q_{ij}}{{\left( {1- t} \right)}^{4- j}}{t^j}} }}{{{{\left[{{\alpha _i}\left( {1-t} \right) + \beta t} \right]}^2}}}。 $

(7)

其中：${Q_{i0}} = {\alpha _i}^2{d_i}, {Q_{i1}} = 2{\alpha _i}^2\left ({3{\Delta _i}-{d_i}} \right), {Q_{i2}} = 3{\alpha _i}{\beta _i}\left ({4{\Delta _i}-{d_i}-{d_{i + 1}}} \right), {\rm{ }}$${Q_{i3}} = 2{\beta _i}^2\left ({3{\Delta _i}-{d_{i + 1}}} \right), {Q_{i4}} = {\beta _i}^2{d_{i + 1}} $。

因为 (7) 式给出的插值函数s(x) 的导数s′(x) 的分母恒为正，所以只需s′(x) 的分子大于0即可。当d_i≥0时，有Q_i1≥0, Q_i4≥0，因此插值函数s(x) 单调的充要条件为：Q_ij≥0, j=1, 2, 3。即

$ \begin{array}{l} {Q_{i1}} = 2{\alpha _i}^2\left( {3{\Delta _i}-{d_i}} \right) \ge 0, {\rm{ }}\\ {Q_{i2}} = 3{\alpha _i}{\beta _i}\left( {4{\Delta _i}-{d_i}-{d_{i + 1}}} \right) \ge 0, {\rm{ }}\\ {\rm{ }}{Q_{i3}} = 2{\beta _i}^2\left( {3{\Delta _i} - {d_{i + 1}}} \right) > 0. \end{array} $

(8)

由 (8) 式得到函数s(x) 单调的充要条件为：

$ {d_i} \le 3{\Delta _i}, {d_{i + 1}} \le 3{\Delta _i}, {d_i} + {d_{i + 1}} \le 4{\Delta _i}。 $

(9)

由文献[8]得出下列结论：

已知单调递增数据f₁≤f₂≤…≤f_n, 导数值d_i≥0, 只要d_i满足 (9) 式时，就存在含有正参数α_i，β_i的有理四次样条插值函数s(x)∈C¹[a, b]，并且是单调递增的。

3 有理四次样条插值函数二阶连续

对 (1) 式求二阶导数并化简可得：

$ s''\left( x \right) = S''\left( t \right){h_i}^{- 2} = \frac{{\sum\limits_{j = 0}^4 {{B_{ij}}{{\left( {1- t} \right)}^{4- j}}{t^j}} }}{{{h_i}{{\left[{{\alpha _i}\left( {1-t} \right) + {\beta _i}t} \right]}^3}}}。 $

(10)

其中：

$ \begin{array}{l} {B_{i0}} = 2{\alpha _i}^3\left( {3{\Delta _i}-2{d_i}} \right)-2{\alpha _i}^2{\beta _i}{d_i}, {\rm{ }}\\ {B_{i1}} = 6{\alpha _i}^2{\beta _i}\left( {3{\Delta _i}-{d_{i + 1}}} \right) + 2{\alpha _i}^3\left( {{d_i} - 3{\Delta _i}} \right) - 8{\alpha _i}^2{\beta _i}{d_i}, {\rm{ }}\\ {B_{i2}} = 6{\alpha _i}^2{\beta _i}\left( {{d_i} - 3{\Delta _i}} \right) + 6{\alpha _i}{\beta _i}^2\left( {3{\Delta _i} - {d_{i + 1}}} \right), \\ {B_{i3}} = 6{\alpha _i}{\beta _i}^2\left( {{d_i} - 3{\Delta _i}} \right) + 2{\beta _i}^3\left( {3{\Delta _i} - {d_{i + 1}}} \right) + 8{\alpha _i}{\beta _i}^2{d_{i + 1}}, {\rm{ }}\\ {B_{i4}} = 2{\beta _i}^3\left( {2{d_{i + 1}} - 3{\Delta _i}} \right) + 2{\alpha _i}{\beta _i}^2{d_{i + 1}}。 \end{array} $

(11)

因此有

$ s''\left( {{x_{i-}}} \right) = \frac{{4\beta _{i-1}^3{d_i} + 2{\alpha _{i-1}}\beta _{i - 1}^2{d_i} - 6\beta _{i - 1}^3{\Delta _{i - 1}}}}{{{h_{i - 1}}\beta _{i - 1}^3}} = \frac{{4{\beta _{i - 1}}{d_i} + 2{\alpha _{i - 1}}{d_i} - 6{\beta _{i - 1}}{\Delta _{i - 1}}}}{{{h_{i - 1}}{\beta _{i - 1}}}}。 $

(12)

$ s''\left( {{x_{i + }}} \right) = \frac{{6\alpha _i^3{\Delta _i}-4\alpha _i^3{d_i}-2\alpha _i^2{\beta _i}{d_i}}}{{{h_i}\alpha _i^3}} = \frac{{6{\alpha _i}{\Delta _i}-4{\alpha _i}{d_i} - 2{\beta _i}{d_i}}}{{{h_i}{\alpha _i}}}。 $

(13)

由函数s(x) 在x_i处二阶连续，可得

$ {d_i} = \frac{{{h_i}{\Delta _{i - 1}} + {h_{i - 1}}{\Delta _i}}}{{\frac{2}{3}\left( {{h_{i - 1}} + {h_i}} \right) + \frac{1}{3}\left( {\frac{{{\beta _i}}}{{{\alpha _i}}}{h_{i - 1}} + \frac{{{\alpha _{i - 1}}}}{{{\beta _{i - 1}}}}{h_i}} \right)}}。 $

(14)

由此得出如下结论：

当插值导数d_i满足 (14) 式时，含有形状控制参数α_i, β_i的有理四次样条插值函数s(x) 在区间[a, b]上二阶连续。^[9]

4 插值曲线的应变能

插值曲线的研究主要应用在物体外形设计上，因此插值曲线不仅要求过插值节点，还要求曲线能够保光顺性、保凹凸性等。研究曲线光顺性会发现光顺准则是最基本的问题，人们在不同的实际问题中，对曲线光顺的要求不同，使得对曲线光顺的认识也不同。一般认为应变能最小的曲线是光顺的，在曲线光顺的研究中常常以曲线能量较小作为约束条件。^[10-12]

插值函数s(x) 在[x₀, x_n]区间上C²-连续曲线的应变能定义如下：

$ \smallint _{{x_0}}^{{x_n}}{\left[{s''\left( x \right)} \right]^2}dx = \sum\limits_{i = 0}^{n - 1} {\smallint _{{x_i}}^{{x_{i + 1}}}{{\left[{s''\left( x \right)} \right]}^2}dx}, $

(15)

积分化简得

$ \sum\limits_{i = 0}^{n- 1} {\smallint _{{x_i}}^{{x_{i + 1}}}{{\left[{s''\left( x \right)} \right]}^2}} dx = \sum\limits_{i = 0}^{n -1} {\frac{{4\left( {{A_i}\alpha _i^2 + {B_i}\beta _i^2 + {C_i}{\alpha _i}{\beta _i}} \right)}}{{5h_i^3{\alpha _i}{\beta _i}}}}, $

(16)

其中：

$ \begin{array}{*{20}{l}} {{A_i} = h_i^2d_{i + 1}^2,}\\ {{B_i} = h_i^2d_i^2,}\\ \begin{array}{l} {C_i} = 15{h_i}{f_i}{d_i} + 15{h_i}{f_i}{d_{i + 1}} + 5h_i^2d_i^2 + 5h_i^2d_{i + 1}^2 - 15{h_i}{f_{i + 1}}{d_i} - \\ 15{h_i}{f_{i + 1}}{d_{i + 1}} + 3h_i^2{d_i}{d_{i + 1}} + 15f_i^2 + 15f_{i + 1}^2 - 30{f_i}{f_{i + 1}}. \end{array} \end{array} $

5 应变能最小的保形有理四次样条插值曲线

为构造光顺的保形有理四次样条插值曲线，以形状控制参数α_i，β_i和节点处的导数d_i为决策变量，以插值曲线应变能最小为目标函数，以插值函数保形以及形状控制参数α_i，β_i和节点处的导数d_i大于0作为约束条件，建立优化模型

$ \smallint _{{x_0}}^{{x_n}}{\left[{s''\left( x \right)} \right]^2}dx = {\rm{min}} $

(17)

$ {\rm{ }}s.t.\;\;\;{\alpha _i}, {\beta _i} > 0, i = 0, 1, \cdots, n-1, $

(18)

$ {\rm{ }}{d_i} > 0i = 0, 1, \cdots, n, $

(19)

$ {d_i} \le 3{\Delta _i}, {d_{i + 1}} \le 3{\Delta _i}, {d_i} + {d_{i + 1}} < 4{\Delta _i}, i = 0, 1, \cdots, n-1. $

(20)

求解此优化模型得最优参数α_i, β_i和d_i，进一步得到光顺的保形有理四次样条插值曲线。

由文献[13-14]知，保形有理四次样条插值函数的误差有以下结论:

定理1 设f(x)∈C²[a, b]，s(x) 是f(x) 如 (1) 式所定义的有理样条插值函数，对给定的α_i, β_i，则x∈[x_i, x_i+1], i=0, 1, …, n-1，有

$ \left| {R[f]} \right| = \left\| {f''(x) -s(x)} \right\| \le \frac{{10}}{9}\left\| {f''(x)} \right\|h_i^2. $

6 数值例子

给出一组单调递增的数据x₁=0, x₂=6, x₃=10, x₄=29.5, x₅=30, f(x₁)=0.01, f(x₂)=15, f(x₃)=15, f(x₄)=25, f(x₅)=30。

由本文方法建立模型求解最优形状参数α_i，β_i和插值节点处导数d_i的值如表 1所示。

依表 1数据绘制出应变能最小且保形的有理四次样条插值曲线 (如图 1)，并与文献[6]中插值法所绘制的插值曲线 (如图 2) 作对比，验证了新方法的有效性。