CN101482979A - 一种光顺优化的nurbs空间曲线曲率连续拼接的cad方法 - Google Patents

Abstract

一种光顺优化的NURBS空间曲线曲率连续拼接的CAD方法，属于空间参数曲线几何造型的CAD领域，其特征在于，初始化阶段输入不连续的两条NURBS曲线C和

；对其中的曲线C通过实行三次延伸操作的方式，使得C与曲线

实现G⁰连续。再调整C曲线延伸部分的控制顶点和相应权重，通过一维黄金分割搜索最小化二阶光顺能量，实现两条曲线的G²连续并且具有最小的光顺能量值。该方法在不改变曲线原有部分的情况下，填补了两条NURBS曲线间的缝隙，并且保证了曲线延伸部分的光顺性最优。

Description

一种光顺优化的NURBS空间曲线曲率连续拼接的CAD方法

技术领域

本发明涉及一种光顺优化的NURBS空间曲线曲率连续拼接的CAD方法。

技术背景

CAD方法在现代制造业中发挥着重要的作用。它可以准确表达较为复杂的自由型曲线曲面，同时能显著改善设计的效率。随着国际标准化组织(ISO)于1991年正式确定把NURBS方法作为定义产品形状的唯一数学方法，越来越多的CAD/CAM系统采用NURBS作为其模型，NURBS在工业产品的外形设计中获得广泛的应用。

工业产品外形设计中的一个重要方面是曲线设计，设计人员通常先设计出最能表现产品外形特点的若干条NURBS曲线，其中的缝隙则希望使用CAD方法进行填补，以满足产品外形的美观要求。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是，提供一种光顺优化的NURBS空间曲线曲率连续拼接的CAD方法，使得在保持用户给定的两条NURBS空间曲线原始部分不变且不引入第三条曲线的前提下，通过三次延伸其中的一条曲线并对此延伸部分进行适当修改，实现两条曲线的曲率连续拼接和拼接部分光顺性最优。

(二)技术方案

本发明公开了一种光顺优化的NURBS空间曲线曲率连续拼接的CAD方法，该方法在计算机中实现，包括如下步骤：

S1：输入不连续的两条NURBS空间曲线C(u)和C(u)，对于曲线C(u)其形式为：

$C\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{n-1}{N}_{i,U}\left(u\right)w_i{P}_i}{{\mathrm{\Σ}}_{i=0}^{n-1}{N}_{i,U}\left(u\right)w_i}$

其中u是曲线C(u)的参数，n是曲线C(u)所含的控制顶点个数，P_i是曲线C(u)中序号为i的控制顶点的坐标值，由用户输入，i＝0，1，...，n-1，采用三维直角坐标表示，w_i是曲线C(u)中序号为i的控制顶点的权重，由用户输入，i＝0，1，...，n-1，N_i，U是定义在节点向量序列U之上的序号为i的B样条基函数，i＝0，1，...，n-1，U是曲线C(u)的节点向量序列：

U＝{u₀，u₁，...，u_n+p}

其中p是曲线C(u)的幂，由用户输入，u₀，u₁，...，u_n+p等为节点向量序列中的节点，其中u₀＝u₁＝…＝u_p＝0，u_n＝u_n+1＝…＝u_n+p＝1，u_p+1，u_p+2，...，u_n-1由用户输入，上述的B样条基函数采用Matlab系统样条工具库中B样条基函数的定义方式，对于曲线C(u)其形式为：

$\stackrel{\‾}{C}\left(\stackrel{\‾}{u}\right)=\frac{{\mathrm{\Σ}}_{i=0}^{\stackrel{\‾}{n}-1}{N}_{i,\stackrel{\‾}{U}}\left(\stackrel{\‾}{u}\right){\stackrel{\‾}{w}}_i{\stackrel{\‾}{P}}_i}{{\mathrm{\Σ}}_{i=0}^{\stackrel{\‾}{n}-1}{N}_{i,\stackrel{\‾}{U}}\left(\stackrel{\‾}{u}\right){\stackrel{\‾}{w}}_i}$

其中u是曲线C(u)的参数，n是曲线C(u)所含的控制顶点个数，P_i是曲线C(u)中序号为i的控制顶点的坐标值，由用户输入，i＝0，1，...，n-1，采用三维直角坐标表示，w_i是曲线C(u)中序号为i的控制顶点的权重，由用户输入，i＝0，1，...，n-1，

是定义在节点向量序列U之上的序号为i的B样条基函数，i＝0，1，...，n-1，U是曲线C(u)的节点向量序列：

$\stackrel{\‾}{U}=\{{\stackrel{\‾}{u}}_0,{\stackrel{\‾}{u}}_1,...,{\stackrel{\‾}{u}}_{\stackrel{\‾}{n}+\stackrel{\‾}{p}}\}$

其中p是曲线C(u)的幂，由用户输入，

等为节点向量序列中的节点，其中 ${\stackrel{\‾}{u}}_0={\stackrel{\‾}{u}}_1=\·\·\·={\stackrel{\‾}{u}}_{\stackrel{\‾}{p}}=0,$ ${\stackrel{\‾}{u}}_{\stackrel{\‾}{n}}={\stackrel{\‾}{u}}_{\stackrel{\‾}{n}+1}=\·\·\·={\stackrel{\‾}{u}}_{\stackrel{\‾}{n}+\stackrel{\‾}{p}}=1,$

由用户输入，上述的B样条基函数采用Matlab系统样条工具库中B样条基函数的定义方式；

S2：对曲线C(u)进行向曲线C(u)方向的延伸，其步骤如下：

S21：计算齐次空间中曲线C(u)的B样条表示形式C^w(u)，计算方法如下：

$C^w\left(u\right)=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{N}_{i,U}\left(u\right)P_i^w$

其中是齐次空间中点P_i的表示形式，i＝0，1，...，n-1，采用四维直角坐标表示，其计算方法如下：

$P_i^w=(w_i{P}_i,w_i)$

计算齐次空间中曲线C(u)的B样条表示形式C^w(u)，计算方法如下：

${\stackrel{\‾}{C}}^w\left(\stackrel{\‾}{u}\right)=\underset{i=0}{\overset{\stackrel{\‾}{n}-1}{\mathrm{\Σ}}}{N}_{i,\stackrel{\‾}{U}}\left(u\right){\stackrel{\‾}{P}}_i^w$

其中是齐次空间中点P_i的表示形式，i＝0，1，...，n-1，采用四维直角坐标表示，其计算方法如下：

${\stackrel{\‾}{P}}_i^w=({\stackrel{\‾}{w}}_i{\stackrel{\‾}{P}}_i,{\stackrel{\‾}{w}}_i),$

S22：在曲线C^w(u)和曲线C^w(u)之间设定三个齐次空间点q₁，q₂和q₃，分别作为三次曲线延伸的目标点，其中点q₁和q₂由用户输入，q₃设定为曲线C^w(u)的第一个控制顶点，即 $q_3={\stackrel{\‾}{P}}_0^w,$

S23：把曲线C^w(u)延伸到所述的点q₁，并且设延伸后的曲线为

，其表示形式为：

$C_1^w\left(u\right)=\underset{i=0}{\overset{n_1-1}{\mathrm{\Σ}}}{N}_{i,U_1}\left(u\right)P_{1,i}^w$

其中n₁是曲线

中所含的控制顶点个数，且有n₁＝n+1，

是曲线

中序号为i的控制顶点的函数值，按照步骤S24所述的方法计算，i＝0，1，...，n₁-1，

是定义在节点向量序列U₁之上的序号为i的B样条基函数，i＝0，1，...，n₁-1，U₁是曲线

的归一化形式的节点向量序列：

其中 $a=1+\frac{\left|\right|P_{1,n_1-1}^w-q_1\left|\right|}{{\mathrm{\Σ}}_{i=p}^{n-1}\left|\right|C^w\left(u_{i+1}\right)-C^w\left(u_i\right)\left|\right|},$ ‖·‖表示齐次空间的欧氏距离，

S24：按下述步骤计算曲线

中序号为i的控制顶点的值

i＝0，1，...，n₁-1，

S241：设置初值 ${\tilde{P}}_j^{-1}=P_{1,j}^w,j=n_1-p,...,n_1+p$

S242：递推计算齐次空间点

$\left\{\begin{array}{c}{\tilde{P}}_j^i={\tilde{P}}_j^{i-1},j=n_1-p,...,n_1-i-2\\ {\tilde{P}}_j^i=\frac{{\tilde{P}}_j^{i-1}-(1-{\mathrm{\α}}_{i,j}){\tilde{P}}_{j-1}^i}{{\mathrm{\α}}_{i,j}},j=n_1-i-1,...,n_1-1\end{array}\right.$

其中 ${\mathrm{\α}}_{i,j}=\frac{u_{n_1}-u_j}{u_{i+j+2}-u_j},$ i＝0，1，...，p-2

S243：计算的最终结果

$\left\{\begin{array}{c}{P}_{1,i}^w=P_i^w,i=0,...,n_1-p-1\\ P_{1,i}^w={\tilde{P}}_i^{p-2},i=n_1-p,...,n_1-1\end{array}\right.$

S25：按照步骤S23、S24所述的方法将所述曲线

延伸到点q₂，并且设延伸后的曲线为

S26：再按照步骤S23、S24所述的方法将所述曲线

延伸到点q₃，并且设延伸后的曲线为

，其表示形式为：

$C_3^w\left(u\right)=\underset{i=0}{\overset{n_3-1}{\mathrm{\Σ}}}{N}_{i,U_3}\left(u\right)P_{3,i}^w$

其中n₃是曲线

所含的控制顶点个数且有n₃＝n+3，是定义在节点向量序列U₃之上的序号为i的B样条基函数，i＝0，1，...，n₃-1，U₃是曲线

的归一化形式的节点向量序列：

$U_3=\{v_0,v_1,...,v_{n_3+p}\},$

其中v₀＝v₁＝…＝v_p＝0， $v_{n_3}=v_{n_3+1}=\·\·\·=v_{n_3+p}=1;$

S27：计算齐次空间曲线在三维空间中的表示形式C₃(u)：

$C_3\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{n_3-1}{N}_{i,U_3}\left(u\right)w_{3,i}{P}_{3,i}}{{\mathrm{\Σ}}_{i=0}^{n_3-1}{N}_{i,U_3}\left(u\right)w_{3,i}}$

其中w_3，i为

的最后一维的数值，i＝0，1，...，n₃-1，设

为

去掉最后一维后的向量，则 $P_{3,i}=\frac{P_{3,i}^{\′}}{w_{3,i}},$ i＝0，1，...，n₃-1；

S3：修改曲线C₃(u)的两个控制顶点P_3，n和P_3，n+1及其相应权重w_3，n和w_3，n+1，使得修改后的曲线C₃(u)与曲线C(u)在其交接处即u＝1和u＝0处实现G²连续，P_3，n和P_3，n+1以及w_3，n和w_3，n+1的计算方法如下：

$w_{3,n+1}={\stackrel{\‾}{w}}_0-\mathrm{\β}\frac{\stackrel{\‾}{\mathrm{\λ}}}{\mathrm{\λ}}({\stackrel{\‾}{w}}_1-{\stackrel{\‾}{w}}_0)$

$P_{3,n+1}=\frac{{\stackrel{\‾}{w}}_0}{w_{3,n+1}}{\stackrel{\‾}{P}}_0-\mathrm{\β}\frac{\stackrel{\‾}{\mathrm{\λ}}}{\mathrm{\λ}{w}_{3,n+1}}({\stackrel{\‾}{w}}_1{\stackrel{\‾}{P}}_1-{\stackrel{\‾}{w}}_0{\stackrel{\‾}{P}}_0)$

$w_{3,n}=-\frac{{\mathrm{\χ}}_{n+2}}{{\mathrm{\χ}}_n}{\stackrel{\‾}{w}}_0+\frac{{\mathrm{\χ}}_{n+1}}{{\mathrm{\χ}}_n}{w}_{3,n+1}+\frac{{\mathrm{\β}}^2}{{\mathrm{\χ}}_n}({\stackrel{\‾}{\mathrm{\χ}}}_2{\stackrel{\‾}{w}}_2-{\stackrel{\‾}{\mathrm{\χ}}}_1{\stackrel{\‾}{w}}_1+{\stackrel{\‾}{\mathrm{\χ}}}_0{\stackrel{\‾}{w}}_0)$

$P_{3,n}=-\frac{{\mathrm{\χ}}_{n+2}{\stackrel{\‾}{w}}_0}{{\mathrm{\χ}}_n{w}_n}{\stackrel{\‾}{P}}_0+\frac{{\mathrm{\χ}}_{n+1}{w}_{n+1}}{{\mathrm{\χ}}_n{w}_n}{P}_{3,n+1}+\frac{{\mathrm{\β}}^2}{{\mathrm{\χ}}_n{w}_n}({\stackrel{\‾}{\mathrm{\χ}}}_2{\stackrel{\‾}{w}}_2{\stackrel{\‾}{P}}_2-{\stackrel{\‾}{\mathrm{\χ}}}_1{\stackrel{\‾}{w}}_1{\stackrel{\‾}{P}}_1+{\stackrel{\‾}{\mathrm{\χ}}}_0{\stackrel{\‾}{w}}_0{\stackrel{\‾}{P}}_0)$

其中β是由步骤S4确定的正实数值，λ，λ，χ_n，χ_n+1，χ_n+2，χ₀，χ₁，χ₂为常数，定义如下：

$\mathrm{\λ}=\frac{p}{1-v_{n_3-1}}$

$\stackrel{\‾}{\mathrm{\λ}}=\frac{\stackrel{\‾}{p}}{{\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+1}}$

${\mathrm{\χ}}_n=\frac{p(p-1)}{(1-v_{n_3-1})(1-v_{n_3-2})}$

${\mathrm{\χ}}_{n+2}=\frac{p(p-1)}{{(1-v_{n_3-1})}^2}$

χ_n+1＝χ_n+χ_n+2

${\stackrel{\‾}{\mathrm{\χ}}}_2=\frac{\stackrel{\‾}{p}(\stackrel{\‾}{p}-1)}{{\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+1}{\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+2}}$

${\stackrel{\‾}{\mathrm{\χ}}}_0=\frac{\stackrel{\‾}{p}(\stackrel{\‾}{p}-1)}{{\left({\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+1}\right)}^2}$

χ₁＝χ₀+χ₂；

S4：确定β的取值使得延伸曲线二阶光顺能量值最小，所述的二阶光顺能量值定义为：

$E=\∫{\left|\right|C_3^{\′\′}\left(u\right)\left|\right|}^2\mathrm{du}$

其中表示C₃(u)的二阶导数，由步骤S3知，C₃(u)的分子分母的各项中只有

和

为β的多项式函数，其余各项均与β无关，故C₃(u)具有如下形式：

$C_3\left(u\right)=\frac{{\mathrm{\β}}^4{v}_1\left(u\right)+{\mathrm{\β}}^2{v}_2\left(u\right)+v_3\left(u\right)}{{\mathrm{\β}}^2{c}_1\left(u\right)+\mathrm{\β}{c}_2\left(u\right)+c_3\left(u\right)}$

其中c₁(u)，c₂(u)，c₃(u)为已知的关于u的分段多项式函数，v₁(u)，v₂(u)，v₃(u)为已知的关于u的分段多项式向量值函数，对C₃(u)关于u求一阶导数可得：

$C_3^{\′}\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^6{\mathrm{\β}}^i{v}_{4+i}\left(u\right)}{{\mathrm{\Σ}}_{j=0}^4{\mathrm{\β}}^j{c}_{4+j}\left(u\right)}$

其中c₄(u)，c₅(u)，...，c₈(u)为已知的关于u的分段多项式函数，v₄(u)，v₅(u)，...，v₁₀(u)为已知的关于u的分段多项式向量值函数，再对

关于u求一阶导数可得：

$C_3^{\′\′}\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{10}{\mathrm{\β}}^i{v}_{11+i}\left(u\right)}{{\mathrm{\Σ}}_{j=0}^8{\mathrm{\β}}^j{c}_{9+j}\left(u\right)}$

其中c₉(u)，c₁₀(u)，...，c₁₇(u)为已知的关于u的分段多项式函数，v₁₁(u)，v₁₂(u)，...，v₂₁(u)为已知的关于u的分段多项式向量值函数，因而

的表达式为：

${\left|\right|C_3^{\′\′}\left(u\right)\left|\right|}^2=C_3^{\′\′}{\left(u\right)}^T{C}_3^{\′\′}\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{20}{\mathrm{\β}}^i{c}_{18+i}\left(u\right)}{{\mathrm{\Σ}}_{j=0}^{16}{\mathrm{\β}}^j{c}_{39+j}\left(u\right)}$

其中c₁₈(u)，c₁₉(u)，...，c₅₅(u)为已知的关于u的分段多项式函数，为分段分式有理函数，因此对

的积分可以得到显式表达式，最后用一维黄金分割搜索法求出最小能量值对应的β。

(三)有益效果

在实际应用中，此方法简单可行。在保持用户给定的两条NURBS空间曲线原始部分不变且不引入第三条曲线的前提下，通过对其中的一条NURBS曲线实行三次延伸操作并对此延伸部分进行适当修改，实现了两条NURBS曲线的曲率连续拼接和光顺性优化。

附图说明

图1为输入的两条原始NURBS曲线C和C；

图2是对曲线C相继实行三次延伸操作后的曲线C₃；

图3为对曲线C₃延伸部分进行修改后与曲线C达到二阶几何连续的结果；

图4为对曲线C₃延伸部分进一步修改后实现延伸部分二阶光顺能量值最小的结果；

图5是整个方法的步骤流程图。

具体实施方式

下面结合附图和具体实施方式对本方法作进一步的详细说明：

本CAD方法以Windows XP中的Visual Studio 2008为开发平台，其具体实施方式如下：

图1是输入的两条原始NURBS曲线C和C。对于曲线C(u)其形式为：

$C\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{n-1}{N}_{i,U}\left(u\right)w_i{P}_i}{{\mathrm{\Σ}}_{i=0}^{n-1}{N}_{i,U}\left(u\right)w_i}$

其中u是曲线C(u)的参数，n是曲线C(u)所含的控制顶点个数，P_i是曲线C(u)中序号为i的控制顶点的坐标值，由用户输入，i＝0，1，...，n-1，采用三维直角坐标表示，w_i是曲线C(u)中序号为i的控制顶点的权重，由用户输入，i＝0，1，...，n-1，N_i，U是定义在节点向量序列U之上的序号为i的B样条基函数，i＝0，1，...，n-1，U是曲线C(u)的节点向量序列：

U＝{u₀，u₁，...，u_n+p}

其中p是曲线C(u)的幂，由用户输入，u₀，u₁，...，u_n+p等为节点向量序列中的节点，其中u₀＝u₁＝…＝u_p＝0，u_n＝u_n+1＝…＝u_n+p＝1，u_p+1，u_p+2，...，u_n-1由用户输入，上述的B样条基函数采用Matlab系统样条工具库中B样条基函数的定义方式，对于曲线C(u)其形式为：

$\stackrel{\‾}{C}\left(\stackrel{\‾}{u}\right)=\frac{{\mathrm{\Σ}}_{i=0}^{\stackrel{\‾}{n}-1}{N}_{i,\stackrel{\‾}{U}}\left(\stackrel{\‾}{u}\right){\stackrel{\‾}{w}}_i{\stackrel{\‾}{P}}_i}{{\mathrm{\Σ}}_{i=0}^{\stackrel{\‾}{n}-1}{N}_{i,\stackrel{\‾}{U}}\left(\stackrel{\‾}{u}\right){\stackrel{\‾}{w}}_i}$

其中u是曲线C(u)的参数，n是曲线C(u)所含的控制顶点个数，P_i是曲线C(u)中序号为i的控制顶点的坐标值，由用户输入，i＝0，1，...，n-1，采用三维直角坐标表示，w_i是曲线C(u)中序号为i的控制顶点的权重，由用户输入，i＝0，1，...，n-1，

是定义在节点向量序列U之上的序号为i的B样条基函数，i＝0，1，...，n-1，U是曲线C(u)的节点向量序列：

$\stackrel{\‾}{U}=\{{\stackrel{\‾}{u}}_0,{\stackrel{\‾}{u}}_1,...,{\stackrel{\‾}{u}}_{\stackrel{\‾}{n}+\stackrel{\‾}{p}}\}$

其中p是曲线C(u)的幂，由用户输入，

等为节点向量序列中的节点，其中 ${\stackrel{\‾}{u}}_0={\stackrel{\‾}{u}}_1=\·\·\·={\stackrel{\‾}{u}}_{\stackrel{\‾}{p}}=0,$ ${\stackrel{\‾}{u}}_{\stackrel{\‾}{n}}={\stackrel{\‾}{u}}_{\stackrel{\‾}{n}+1}=\·\·\·={\stackrel{\‾}{u}}_{\stackrel{\‾}{n}+\stackrel{\‾}{p}}=1,$

由用户输入，上述的B样条基函数采用Matlab系统样条工具库中B样条基函数的定义方式。

图2是对曲线C相继实行三次延伸操作后的结果，其具体方法如下：

(1)计算齐次空间中曲线C(u)的B样条表示形式C^w(u)，计算方法如下：

$C^w\left(u\right)=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{N}_{i,U}\left(u\right)P_i^w$

其中

是齐次空间中点P_i的表示形式，i＝0，1，...，n-1，采用四维直角坐标表示，其计算方法如下：

$P_i^w=(w_i{P}_i,w_i)$

计算齐次空间中曲线C(u)的B样条表示形式C^w(u)，计算方法如下：

${\stackrel{\‾}{C}}^w\left(\stackrel{\‾}{u}\right)=\underset{i=0}{\overset{\stackrel{\‾}{n}-1}{\mathrm{\Σ}}}{N}_{i,\stackrel{\‾}{U}}\left(u\right){\stackrel{\‾}{P}}_i^w$

其中

是齐次空间中点Pi的表示形式，i＝0，1，...，n-1，采用四维直角坐标表示，其计算方法如下：

${\stackrel{\‾}{P}}_i^w=({\stackrel{\‾}{w}}_i{\stackrel{\‾}{P}}_i,{\stackrel{\‾}{w}}_i),$

(2)在曲线C^w(u)和曲线C^w(u)之间设定三个齐次空间点q₁，q₂和q₃，分别作为三次曲线延伸的目标点，其中点q₁和q₂由用户输入，q₃设定为曲线C^w的第一个控制顶点，即 $q_3={\stackrel{\‾}{P}}_0^w,$

(3)把曲线C^w(u)延伸到所述的点q₁，并且设延伸后的曲线为

其表示形式为：

$C_1^w\left(u\right)=\underset{i=0}{\overset{n_1-1}{\mathrm{\Σ}}}{N}_{i,U_1}\left(u\right)P_{1,i}^w$

其中n₁是曲线

中所含的控制顶点个数，且有n₁＝n+1，

是定义在节点向量序列U₁之上的序号为i的B样条基函数，i＝0，1，...，n₁-1，U₁是曲线

的归一化形式的节点向量序列：

其中 $a=1+\frac{\left|\right|P_{1,n_1-1}^w-q_1\left|\right|}{{\mathrm{\Σ}}_{i=p}^{n-1}\left|\right|C^w\left(u_{i+1}\right)-C^w\left(u_i\right)\left|\right|},$ ‖·‖表示齐次空间的欧氏距离，

是曲线

中序号为i的控制顶点的函数值，i＝0，1，...，n₁-1，其计算方法如下：

(a)设置初值 ${\tilde{P}}_j^{-1}=P_{1,j}^w,$ j＝n₁-p，...，n₁+p

(b)递推计算齐次空间点

$\left\{\begin{array}{c}{\tilde{P}}_j^i={\tilde{P}}_j^{i-1},j=n_1-p,...,n_1-i-2\\ {\tilde{P}}_j^i=\frac{{\tilde{P}}_j^{i-1}-(1-{\mathrm{\α}}_{i,j}){\tilde{P}}_{j-1}^i}{{\mathrm{\α}}_{i,j}},j=n_1-i-1,...,n_1-1\end{array}\right.$

其中 ${\mathrm{\α}}_{i,j}=\frac{u_{n_1}-u_j}{u_{i+j+2}-u_j},$ i＝0，1，...，p-2

(c)计算的最终结果

$\left\{\begin{array}{c}{P}_{1,i}^w=P_i^w,i=0,...,n_1-p-1\\ P_{1,i}^w={\tilde{P}}_i^{p-2},i=n_1-p,...,n_1-1\end{array}\right.$

(4)按照(3)中所述的方法将所述曲线

延伸到点q₂，并且设延伸后的曲线为

(5)再按照(3)中所述的方法将所述曲线

延伸到点q₃，并且设延伸后的曲线为

，其表示形式为：

$C_3^w\left(u\right)=\underset{i=0}{\overset{n_3-1}{\mathrm{\Σ}}}{N}_{i,U_3}\left(u\right)P_{3,i}^w$

其中n₃是曲线

所含的控制顶点个数且有n₃＝n+3，

是定义在节点向量序列U₃之上的序号为i的B样条基函数，i＝0，1，...，n₃-1，U₃是曲线

的归一化形式的节点向量序列：

$U_3=\{v_0,v_1,...,v_{n_3+p}\},$

其中v₀＝v₁＝…＝v_p＝0， $v_{n_3}=v_{n_3+1}=\·\·\·=v_{n_3+p}=1;$

(6)计算齐次空间曲线

在三维空间中的表示形式C₃(u)：

$C_3\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{n_3-1}{N}_{i,U_3}\left(u\right)w_{3,i}{P}_{3,i}}{{\mathrm{\Σ}}_{i=0}^{n_3-1}{N}_{i,U_3}\left(u\right)w_{3,i}}$

其中w_3，i为

的最后一维的数值，i＝0，1，...，n₃-1，设

为

去掉最后一维后的向量，则 $P_{3,i}=\frac{P_{3,i}^{\′}}{w_{3,i}},$ i＝0，1，...，n₃-1。

图3为对曲线C₃延伸部分进行修改后与曲线C达到曲率连续的结果，其具体方法如下：

修改曲线C₃(u)的两个控制顶点P_3，n和P_3，n+1及其相应权重w_3，n和w_3，n+1，使得修改后的曲线C₃(u)与曲线C(u)在其交接处即u＝1和u＝0处实现G²连续，P_3，n和P_3，n+1以及w_3，n和w_3，n+1的计算方法如下：

$w_{3,n+1}={\stackrel{\‾}{w}}_0-\mathrm{\β}\frac{\stackrel{\‾}{\mathrm{\λ}}}{\mathrm{\λ}}({\stackrel{\‾}{w}}_1-{\stackrel{\‾}{w}}_0)$

$P_{3,n+1}=\frac{{\stackrel{\‾}{w}}_0}{w_{3,n+1}}{\stackrel{\‾}{P}}_0-\mathrm{\β}\frac{\stackrel{\‾}{\mathrm{\λ}}}{\mathrm{\λ}{w}_{3,n+1}}({\stackrel{\‾}{w}}_1{\stackrel{\‾}{P}}_1-{\stackrel{\‾}{w}}_0{\stackrel{\‾}{P}}_0)$

$w_{3,n}=-\frac{{\mathrm{\χ}}_{n+2}}{{\mathrm{\χ}}_n}{\stackrel{\‾}{w}}_0+\frac{{\mathrm{\χ}}_{n+1}}{{\mathrm{\χ}}_n}{w}_{3,n+1}+\frac{{\mathrm{\β}}^2}{{\mathrm{\χ}}_n}({\stackrel{\‾}{\mathrm{\χ}}}_2{\stackrel{\‾}{w}}_2-{\stackrel{\‾}{\mathrm{\χ}}}_1{\stackrel{\‾}{w}}_1+{\stackrel{\‾}{\mathrm{\χ}}}_0{\stackrel{\‾}{w}}_0)$

$P_{3,n}=-\frac{{\mathrm{\χ}}_{n+2}{\stackrel{\‾}{w}}_0}{{\mathrm{\χ}}_n{w}_n}{\stackrel{\‾}{P}}_0+\frac{{\mathrm{\χ}}_{n+1}{w}_{n+1}}{{\mathrm{\χ}}_n{w}_n}{P}_{3,n+1}+\frac{{\mathrm{\β}}^2}{{\mathrm{\χ}}_n{w}_n}({\stackrel{\‾}{\mathrm{\χ}}}_2{\stackrel{\‾}{w}}_2{\stackrel{\‾}{P}}_2-{\stackrel{\‾}{\mathrm{\χ}}}_1{\stackrel{\‾}{w}}_1{\stackrel{\‾}{P}}_1+{\stackrel{\‾}{\mathrm{\χ}}}_0{\stackrel{\‾}{w}}_0{\stackrel{\‾}{P}}_0)$

其中β是待定的正实数值，λ，λ，χ_n，χ_n+1，χ_n+2，χ₀，χ₁，χ₂为常数，定义如下：

$\mathrm{\λ}=\frac{p}{1-v_{n_3-1}}$

$\stackrel{\‾}{\mathrm{\λ}}=\frac{\stackrel{\‾}{p}}{{\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+1}}$

${\mathrm{\χ}}_n=\frac{p(p-1)}{(1-v_{n_3-1})(1-v_{n_3-2})}$

${\mathrm{\χ}}_{n+2}=\frac{p(p-1)}{{(1-v_{n_3-1})}^2}$

χ_n+1＝χ_n+χ_n+2

${\stackrel{\‾}{\mathrm{\χ}}}_2=\frac{\stackrel{\‾}{p}(\stackrel{\‾}{p}-1)}{{\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+1}{\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+2}}$

${\stackrel{\‾}{\mathrm{\χ}}}_0=\frac{\stackrel{\‾}{p}(\stackrel{\‾}{p}-1)}{{\left({\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+1}\right)}^2}$

χ₁＝χ₀+χ₂

图4为对曲线C₃延伸部分进一步修改后实现延伸部分二阶光顺能量值最小的结果，其具体方法如下：

确定β的取值使得延伸曲线二阶光顺能量值最小，所述的二阶光顺能量值定义为：

$E=\∫{\left|\right|C_3^{\′\′}\left(u\right)\left|\right|}^2\mathrm{du}$

其中

表示C₃(u)的二阶导数，

由于C₃(u)的分子分母的各项中只有

和为β的多项式函数，其余各项均与β无关，故C₃(u)具有如下形式：

$C_3\left(u\right)=\frac{{\mathrm{\β}}^4{v}_1\left(u\right)+{\mathrm{\β}}^2{v}_2\left(u\right)+v_3\left(u\right)}{{\mathrm{\β}}^2{c}_1\left(u\right)+\mathrm{\β}{c}_2\left(u\right)+c_3\left(u\right)}$

其中c₁(u)，c₂(u)，c₃(u)为已知的关于u的分段多项式函数，v₁(u)，v₂(u)，v₃(u)为已知的关于u的分段多项式向量值函数，对C₃(u)关于u求一阶导数可得：

$C_3^{\′}\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^6{\mathrm{\β}}^i{v}_{4+i}\left(u\right)}{{\mathrm{\Σ}}_{j=0}^4{\mathrm{\β}}^j{c}_{4+j}\left(u\right)}$

其中c₄(u)，c₅(u)，...，c₈(u)为已知的关于u的分段多项式函数，v₄(u)，v₅(u)，...，v₁₀(u)为已知的关于u的分段多项式向量值函数，再对

关于u求一阶导数可得：

$C_3^{\′\′}\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{10}{\mathrm{\β}}^i{v}_{11+i}\left(u\right)}{{\mathrm{\Σ}}_{j=0}^8{\mathrm{\β}}^j{c}_{9+j}\left(u\right)}$

其中c₉(u)，c₁₀(u)，...，c₁₇(u)为已知的关于u的分段多项式函数，v₁₁(u)，v₁₂(u)，...，v₂₁(u)为已知的关于u的分段多项式向量值函数，因而

的表达式为：

${\left|\right|C_3^{\′\′}\left(u\right)\left|\right|}^2=C_3^{\′\′}{\left(u\right)}^T{C}_3^{\′\′}\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{20}{\mathrm{\β}}^i{c}_{18+i}\left(u\right)}{{\mathrm{\Σ}}_{j=0}^{16}{\mathrm{\β}}^j{c}_{39+j}\left(u\right)}$

其中c₁₈(u)，c₁₉(u)，...，c₅₅(u)为已知的关于u的分段多项式函数，

为分段分式有理函数，使用Mathematica软件可以求出二阶光顺能量值E关于β的表达式，并最终求出最小能量值对应的β。

图5是整个方法的步骤流程图，初始化阶段输入不连续的两条NURBS曲线C和C；对其中的曲线C通过实行三次延伸操作的方式，使得C与曲线C实现G⁰连续。再调整C曲线延伸部分的控制顶点和相应权重，通过一维黄金分割搜索最小化二阶光顺能量，实现两条曲线的G²连续并且具有最小的光顺能量值。

Claims (1)

1.一种光顺优化的NURBS空间曲线曲率连续拼接的CAD方法，其特征在于，所述方法是在计算机上依次按照如下步骤实现的：

S1：输入不连续的两条NURBS空间曲线C(u)和C(u)，对于曲线C(u)其形式为：

$C\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{n-1}{N}_{i,U}\left(u\right)w_i{P}_i}{{\mathrm{\Σ}}_{i=0}^{n-1}{N}_{i,U}\left(u\right)w_i}$

其中u是曲线C(u)的参数，n是曲线C(u)所含的控制顶点个数，P_i是曲线C(u)中序号为i的控制顶点的坐标值，由用户输入，i＝0，1，...，n-1，采用三维直角坐标表示，w_i是曲线C(u)中序号为i的控制顶点的权重，由用户输入，i＝0，1，...，n-1，N_i，U是定义在节点向量序列U之上的序号为i的B样条基函数，i＝0，1，...，n-1，U是曲线C(u)的节点向量序列：

U＝{u₀，u₁，...，u_n+p}

其中p是曲线C(u)的幂，由用户输入，u₀，u₁，...，u_n+p等为节点向量序列中的节点，其中u₀＝u₁＝…＝u_p＝0，u_n＝u_n+1＝…＝u_n+p＝1，u_p+1，u_p+2，...，u_n-1由用户输入，上述的B样条基函数采用Matlab系统样条工具库中B样条基函数的定义方式；对于曲线C(u)其形式为：

$\stackrel{\‾}{C}\left(\stackrel{\‾}{u}\right)=\frac{{\mathrm{\Σ}}_{i=0}^{\stackrel{\‾}{n}-1}{N}_{i,\stackrel{\‾}{U}}\left(\stackrel{\‾}{u}\right){\stackrel{\‾}{w}}_i{\stackrel{\‾}{P}}_i}{{\mathrm{\Σ}}_{i=0}^{\stackrel{\‾}{n}-1}{N}_{i,\stackrel{\‾}{U}}\left(\stackrel{\‾}{u}\right){\stackrel{\‾}{w}}_i}$

其中u是曲线C(u)的参数，n是曲线C(u)所含的控制顶点个数，P_i是曲线C(u)中序号为i的控制顶点的坐标值，由用户输入，i＝0，1，...，n-1，采用三维直角坐标表示，w_i是曲线C(u)中序号为i的控制顶点的权重，由用户输入，i＝0，1，...，n-1，

是定义在节点向量序列U之上的序号为i的B样条基函数，i＝0，1，...，n-1，U是曲线C(u)的节点向量序列：

$\stackrel{\‾}{U}=\{{\stackrel{\‾}{u}}_0,{\stackrel{\‾}{u}}_1,...,{\stackrel{\‾}{u}}_{\stackrel{\‾}{n}+\stackrel{\‾}{p}}\}$

其中p是曲线C(u)的幂，由用户输入，

等为节点向量序列中的节点，其中 ${\stackrel{\‾}{u}}_0={\stackrel{\‾}{u}}_1=\·\·\·={\stackrel{\‾}{u}}_{\stackrel{\‾}{p}}=0,$ ${\stackrel{\‾}{u}}_{\stackrel{\‾}{n}}={\stackrel{\‾}{u}}_{\stackrel{\‾}{n}+1}=\·\·\·={\stackrel{\‾}{u}}_{\stackrel{\‾}{n}+\stackrel{\‾}{p}}=1,$ 由用户输入，上述的B样条基函数采用Matlab系统样条工具库中B样条基函数的定义方式；

S2：对曲线C(u)进行向曲线C(u)方向的延伸，其步骤如下：

S21：计算齐次空间中曲线C(u)的B样条表示形式C^w(u)，计算方法如下：

$C^w\left(u\right)=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{N}_{i,U}\left(u\right)P_i^w$

其中

是齐次空间中点P_i的表示形式，i＝0，1，...，n-1，采用四维直角坐标表示，其计算方法如下：

$P_i^w=(w_i{P}_i,w_i)$

计算齐次空间中曲线C(u)的B样条表示形式C^w(u)，计算方法如下：

${\stackrel{\‾}{C}}^w\left(\stackrel{\‾}{u}\right)=\underset{i=0}{\overset{\stackrel{\‾}{n}-1}{\mathrm{\Σ}}}{N}_{i,\stackrel{\‾}{U}}\left(u\right){\stackrel{\‾}{P}}_i^w$

其中

是齐次空间中点P_i的表示形式，i＝0，1，...，n-1，采用四维直角坐标表示，其计算方法如下：

${\stackrel{\‾}{P}}_i^w=({\stackrel{\‾}{w}}_i{\stackrel{\‾}{P}}_i,{\stackrel{\‾}{w}}_i),$

S22：在曲线C^w(u)和曲线C^w(u)之间设定三个齐次空间点q₁，q₂和q₃，分别作为三次曲线延伸的目标点，其中点q₁和q₂由用户输入，q₃设定为曲线C^w(u)的第一个控制顶点，即 $q_3={\stackrel{\‾}{P}}_0^w,$

S23：把曲线C^w(u)延伸到所述的点q₁，并且设延伸后的曲线为

其表示形式为：

$C_1^w\left(u\right)=\underset{i=0}{\overset{n_1-1}{\mathrm{\Σ}}}{N}_{i,U_1}\left(u\right)P_{1,i}^w$

其中n₁是曲线

中所含的控制顶点个数，且有n₁＝n+1，

是曲线

中序号为i的控制顶点的函数值，按照步骤S24所述的方法计算，i＝0，1，...，n₁-1，

是定义在节点向量序列U₁之上的序号为i的B样条基函数，i＝0，1，...，n₁-1，U₁是曲线

的归一化形式的节点向量序列：

其中 $a=1+\frac{\left|\right|P_{1,n_1-1}^w-q_1\left|\right|}{{\mathrm{\Σ}}_{i=p}^{n-1}\left|\right|C^w\left(u_{i+1}\right)-C^w\left(u_i\right)\left|\right|},$ ‖·‖表示齐次空间的欧氏距离，

S24：按下述步骤计算曲线中序号为i的控制顶点的值i＝0，1，...，n₁-1，

S241：设置初值 ${\tilde{P}}_j^{-1}=P_{1,j}^w,$ j＝n₁-p，...，n₁+p

S242：递推计算齐次空间点

$\left\{\begin{array}{c}{\tilde{P}}_j^i={\tilde{P}}_j^{i-1},j=n_1-p,...,n_1-i-2\\ {\tilde{P}}_j^i=\frac{{\tilde{P}}_j^{i-1}-(1-{\mathrm{\α}}_{i,j}){\tilde{P}}_{j-1}^i}{{\mathrm{\α}}_{i,j}},j=n_1-i-1,...,n_1-1\end{array}\right.$

其中 ${\mathrm{\α}}_{i,j}=\frac{u_{n_1}-u_j}{u_{i+j+2}-u_j},$ i＝0，1，...，p-2

S243：计算

的最终结果

$\left\{\begin{array}{c}{P}_{1,i}^w=P_i^w,i=0,...,n_1-p-1\\ P_{1,i}^w={\tilde{P}}_i^{p-2},i=n_1-p,...,n_1-1\end{array}\right.$

S25：按照步骤S23、S24所述的方法将所述曲线

延伸到点q₂，并且设延伸后的曲线为

S26：再按照步骤S23、S24所述的方法将所述曲线

延伸到点q₃，并且设延伸后的曲线为

其表示形式为：

$C_3^w\left(u\right)=\underset{i=0}{\overset{n_3-1}{\mathrm{\Σ}}}{N}_{i,U_3}\left(u\right)P_{3,i}^w$

其中n₃是曲线

所含的控制顶点个数且有n₃＝n+3，是定义在节点向量序列U₃之上的序号为i的B样条基函数，i＝0，1，...，n₃-1，U₃是曲线

的归一化形式的节点向量序列：

$U_3=\{v_0,v_1,...,v_{n_3+p}\},$

其中v₀＝v₁＝…＝v_p＝0， $v_{n_3}=v_{n_3+1}=\·\·\·=v_{n_3+p}=1;$

S27：计算齐次空间曲线

在三维空间中的表示形式C₃(u)：

$C_3\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{n_3-1}{N}_{i,U_3}\left(u\right)w_{3,i}{P}_{3,i}}{{\mathrm{\Σ}}_{i=0}^{n_3-1}{N}_{i,U_3}\left(u\right)w_{3,i}}$

其中w_3，i为

的最后一维的数值，i＝0，1，...，n₃-1，设

为

去掉最后一维后的向量，则 $P_{3,i}=\frac{P_{3,i}^{\′}}{w_{3,i}},$ i＝0，1，...，n₃-1；

S3：修改曲线C₃(u)的两个控制顶点P_3，n和P_3，n+1及其相应权重w_3，n和w_3，n+1，使得修改后的曲线C₃(u)与曲线C(u)在其交接处即u＝1和u＝0处实现G²连续，P_3，n和P_3，n+1以及w_3，n和w_3，n+1的计算方法如下：

$w_{3,n+1}={\stackrel{\‾}{w}}_0-\mathrm{\β}\frac{\stackrel{\‾}{\mathrm{\λ}}}{\mathrm{\λ}}({\stackrel{\‾}{w}}_1-{\stackrel{\‾}{w}}_0)$

$P_{3,n+1}=\frac{{\stackrel{\‾}{w}}_0}{w_{3,n+1}}{\stackrel{\‾}{P}}_0-\mathrm{\β}\frac{\stackrel{\‾}{\mathrm{\λ}}}{{\mathrm{\λw}}_{3,n+1}}({\stackrel{\‾}{w}}_1{\stackrel{\‾}{P}}_1-{\stackrel{\‾}{w}}_0{\stackrel{\‾}{P}}_0)$

$w_{3,n}=-\frac{{\mathrm{\χ}}_{n+2}}{{\mathrm{\χ}}_n}{\stackrel{\‾}{w}}_0+\frac{{\mathrm{\χ}}_{n+1}}{{\mathrm{\χ}}_n}{w}_{3,n+1}+\frac{{\mathrm{\β}}^2}{{\mathrm{\χ}}_n}({\stackrel{\‾}{\mathrm{\χ}}}_2{\stackrel{\‾}{w}}_2-{\stackrel{\‾}{\mathrm{\χ}}}_1{\stackrel{\‾}{w}}_1+{\stackrel{\‾}{\mathrm{\χ}}}_0{\stackrel{\‾}{w}}_0)$

$P_{3,n}=-\frac{{\mathrm{\χ}}_{n+2}{\stackrel{\‾}{w}}_0}{{\mathrm{\χ}}_n{w}_n}{\stackrel{\‾}{P}}_0+\frac{{\mathrm{\χ}}_{n+1}{w}_{n+1}}{{\mathrm{\χ}}_n{w}_n}{P}_{3,n+1}+\frac{{\mathrm{\β}}^2}{{\mathrm{\χ}}_n{w}_n}({\stackrel{\‾}{\mathrm{\χ}}}_2{\stackrel{\‾}{w}}_2{\stackrel{\‾}{P}}_2-{\stackrel{\‾}{\mathrm{\χ}}}_1{\stackrel{\‾}{w}}_1{\stackrel{\‾}{P}}_1+{\stackrel{\‾}{\mathrm{\χ}}}_0{\stackrel{\‾}{w}}_0{\stackrel{\‾}{P}}_0)$

其中β是由步骤S4确定的正实数值，λ，λ，x_n，x_n+1，x_n+2，x₀，x₁，x₂为常数，定义如下：

$\mathrm{\λ}=\frac{p}{1-v_{n_3-1}}$

$\stackrel{\‾}{\mathrm{\λ}}=\frac{\stackrel{\‾}{p}}{{\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+1}}$

${\mathrm{\χ}}_n=\frac{p(p-1)}{(1-v_{n_3-1})(1-v_{n_3-2})}$

${\mathrm{\χ}}_{n+2}=\frac{p(p-1)}{{(1-v_{n_3-1})}^2}$

x_n+1＝x_n+x_n+2

${\stackrel{\‾}{\mathrm{\χ}}}_2=\frac{\stackrel{\‾}{p}(\stackrel{\‾}{p}-1)}{{\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+1}{\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+2}}$

${\stackrel{\‾}{\mathrm{\χ}}}_0=\frac{\stackrel{\‾}{p}(\stackrel{\‾}{p}-1)}{{\left({\stackrel{\‾}{u}}_{\stackrel{\‾}{p}+1}\right)}^2}$

x₁＝x₀+x₂；

S4：确定β的取值使得延伸曲线二阶光顺能量值最小，所述的二阶光顺能量值定义为：

$E=\∫{\left|\right|C_3^{\′\′}\left(u\right)\left|\right|}^2\mathrm{du}$

其中

表示C₃(u)的二阶导数；由步骤S3知，C₃(u)的分子分母的各项中只有

和为β的多项式函数，其余各项均与β无关，故C₃(u)具有如下形式：

$C_3\left(u\right)=\frac{{\mathrm{\β}}^4{v}_1\left(u\right)+{\mathrm{\β}}^2{v}_2\left(u\right)+v_3\left(u\right)}{{\mathrm{\β}}^2{c}_1\left(u\right)+\mathrm{\β}{c}_2\left(u\right)+c_3\left(u\right)}$

其中c₁(u)，c₂(u)，c₃(u)为已知的关于u的分段多项式函数，v₁(u)，v₂(u)，v₃(u)为已知的关于u的分段多项式向量值函数，对C₃(u)关于u求一阶导数可得：

$C_3^{\′}\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^6{\mathrm{\β}}^i{v}_{4+i}\left(u\right)}{{\mathrm{\Σ}}_{j=0}^4{\mathrm{\β}}^j{c}_{4+j}\left(u\right)}$

其中c₄(u)，c₅(u)，...，c₈(u)为已知的关于u的分段多项式函数，v₄(u)，v₅(u)，...，v₁₀(u)为已知的关于u的分段多项式向量值函数，再对

关于u求一阶导数可得：

$C_3^{\′\′}\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{10}{\mathrm{\β}}^i{v}_{11+i}\left(u\right)}{{\mathrm{\Σ}}_{j=0}^8{\mathrm{\β}}^j{c}_{9+j}\left(u\right)}$

其中c₉(u)，c₁₀(u)，...，c₁₇(u)为已知的关于u的分段多项式函数，v₁₁(u)，v₁₂(u)，...，v₂₁(u)为已知的关于u的分段多项式向量值函数，因而的表达式为：

${\left|\right|C_3^{\′\′}\left(u\right)\left|\right|}^2=C_3^{\′\′}{\left(u\right)}^T{C}_3^{\′\′}\left(u\right)=\frac{{\mathrm{\Σ}}_{i=0}^{20}{\mathrm{\β}}^i{c}_{18+i}\left(u\right)}{{\mathrm{\Σ}}_{j=0}^{16}{\mathrm{\β}}^j{c}_{39+j}\left(u\right)}$

其中c₁₈(u)，c₁₉(u)，...，c₅₅(u)为已知的关于u的分段多项式函数，

为分段分式有理函数，因此对

的积分可以得到显式表达式，最后用一维黄金分割搜索法求出最小能量值对应的β。

Images