三次样条插值函数的构造与Matlab实现.doc

《三次样条插值函数的构造与Matlab实现.doc》由会员分享，可在线阅读，更多相关《三次样条插值函数的构造与Matlab实现.doc（3页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、兵工 自动 化 自 动测量 与控 制 O. I. Automation 2006 年第 25 卷 第 11 期 Automatic Measurement a nd Contr ol 2006, Vol. 25, No. 11 文章 编号： 1006-1576 （ 2006） 11-0076-03 三次样条插值函数的构造与 Matlab 实现 许 小勇 ，钟 太 勇 （ 1. 云南 民 族大 学 数 学 与计 算机 科 学学 院，云南 昆 明 650031； 2. 郧 阳 师范 高等 专 科学 校 数 学 系，湖北 丹江 口 442700） 摘要 ： 三 次 样条 插值 函 数边 界条 件 由

2、实 际问 题 对三 次样 条 插值 在端 点 的状 态要 求 给出 。以 第 1 边 界条 件为 例 ， 用 节 点处 二阶 导 数表 示三 次 样条 插值 函 数， 用追 赶 法求 解相 关 方程 组。 通 过 Matlab 编 制三 次 样条 函数 的 通用 程序 ， 可 直 接显 示各 区 间段 三次 样 条函 数体 表 达式 ，计 算 出已 给点 插 值并 显示 各 区间 分段 曲 线图 。 关键 词： 三 次样 条； 插 值函 数； Ma tla b 程序 中图 分类 号 ： O242.1 文 献标 识 码： A Constru ctio n an d Real izatio n o

3、f Cu bic Sp line Interpol atio n Fun ctio n XU Xiao-yong , ZHONG Ta i-yong (1. Sc hool of Ma the matics & Compute r Science , Yunnan Nationalities Univer sity, Kunming, 650031, China; 2. De pt. of Ma the matics, Yunyang Teac he r Colle ge , Danjia ngkou 442700, China) Abstrac t: The boundar y condit

4、ions of cubic spline were give n by the pra ctical state at the extreme points of discre te function. Ta king the f ir st boundary condition a s an example, it use d the se cond order derivative of the junctions to repre sent the c ubic spline inter polation func tion, a nd solve d the r elative equ

5、ations with c ha sing method. The unive rsal procedure of cubic spline was re aliz ed unde r the environment of Ma tlab, which c ould compute the inte rpolation point and dir ectly display the exac t expre ssion of function and segmental c ur ve in each inte rval. Ke ywords: Cubic spline; Interpolat

6、ion function; Matlab pr ogr am 0 引言 式 (2)、 (3)共给出 n 3(n-2) 4n-6 个条 件， 称为自然边界条件 。 1 三次样条函数的定义及特征 第 3 类边界条件是周期 性条件，如果 y f(x) 1,2, ,n-1）上都是 不高于三次的多项式（为了与其 对应 j 从 1 开始，在 Matlab 中元素脚标从 1 开始）。 2 三次样条插值函数的推导过程 j j +1 已知，设 S(xj )= Mj ， S(xj+1)=Mj+1，则 S (x) 的表达 j j j j j j j +1 j 收稿日期 ： 2006-07 -21；修回日期： 20

7、06-09-09 作者简介 ：许小勇（ 1983-），男，江西人，云南民族 大学在读硕士，从事 智能计算、小波分析、 图像处理研究。 76 1 1,2 1 1, 2 需要待定 4(n-1) 个系数，因此要唯 一确定三次插值 分段低次样条插值虽然计算简单、稳定性 好、 函数，还要附加2 个边界条件。通 常由实际问题对 收敛性有保证且易在电子计算机 上实现，但只能保 三次样条插值在端点的状态要求给出 。常用边界的 证各小段曲线在连接处的连续性 ，不能保证整件曲 条件有以下3 类。 线的光滑性。利用样条插值，既 可保持分段低次插 第1 类边界条件：给定端点处 的一阶导数值， 值多项式，又可提高插值函

8、数光 滑性。故给出分段 S(x) y ， S(x ) y 。 三次样条插值的构造过程、算法步骤，利用 Matlab 第2 类边界条件：给定端点处 的二阶导数值， 软件编写三次样条插值函数通用 程序，并通过数值 S(x )y ， S(x ) y 。特 殊情况 y y 0， 算例证明程序的正确性。 是以 b-a 为周期的函数，于是 S(x) 在端点处满足条 定义： 设a,b 上有插值节点， a x x 件 S(x 0) S(x -0)， S(x 0) S(x -0)。 x b，对 应函 数值 为 y ,y , y 。 若函 数 S(x) 满 下以第 1 边界条件为例，利用 节点处二阶导数 足 S(

9、x ) y （ j 1,2, ,n）， S(x) 在 x ,x （ j 来表示三次样条插值函数，给出具体的推 导过程。 当 S(x) 在 a,b 具有二阶连续导 数。则称 S(x) 为三 j j+1 次 样 条 插 值 函 数 。 要 求 S (x) 只 需 在 每 个 子 区 间 x ,x 上 确定 1 个三次多项式， 设为： 如果 S(x) 在 x ,x （ j=1,2, ,n-1）两端点上的值 3 2 S( x) a x b x c x d , (j 1,2, ,n-1) ( 1) x x x x j j j j 式为 ： S (x ) M M ， 其 中 h j j1 S(x) y ,

10、 S( x -0)=S(x 0) , (j=2, ,n-1) (2) h h S(x -0)=S(x 0) , S(x -0) =S(x 0) , (j=2, ,n-1) (3) 兵工 自动 化 自 动测量 与控 制 O. I. Automation 2006 年第 25 卷 第 11 期 Automatic Measurement a nd Contr ol 2006, Vol. 25, No. 11 j j j+1 j +1 n n 1 n n j 方程 组 (11)的系数 矩阵 是三对 角阵 且是对 角占 得的解代入 (7)， 就可以构造 a, b上的插值函 数。 6h j 6h j 6

11、 h j 的表达式不 同，但由于要保证在 结点处的连续性， j j1 j j Sj (xj 0) Sj1(xj 0) 程，源程序如下： disp(x0 为 插值 点 ) 整 理得： x0 for i=2:N-1 77 个任意常数 A ， B 的 S(x) 表达式。对 S(x) 求两次 1 ， 则 M 2M (1- )M d （ j 3 3 (x x) ( x x ) ， h 2,3, ,n-1）。把边界条 件 S(x ) y , S (x ) y 代 积分得 S(x) M M A x B 6h 6h 入 S(x)得： x -x （ j =1,2, ,n-1），根据插值 条件 S (x ) y

12、， 2M M (10) M 2M 3 (x x ) 6 y y M A x B y j j j j j 其中： 6h h h 3 (x x ) j1 j (5) y n n 6h 将方程组写成矩阵 形式， j1 j 由式 (5) 、 式 (4)得 ， j j j j1 j 2 1 M 1 1 2 M d y y M M 。将 A 代入 (5)可得， 则 O O O M M h 6 2 M d 2 2 M ，则： j j1 j j j1 6 h 6 j1 j j1 j j1 优阵，故存在 唯 一解，可用追赶法进 行求解，将求 Ax B h x y h j y y M M (6) h x 3 算法

13、步骤 j 输 入 n 个插值结点， a x x x b， M x x M x x 2 2 y h h y 对应函数值为 y ,y , ,y ，边界条件 y ,y ，待求插 j j 值点 x 。 因此可得： 3 3 计 算 h x -x （ j 1,2, ,n-1）。 (x x) (x x ) M x x 2 S(x)M M y h (7)h y y y y 2 ， (j=2, ,n-1)。 j j j1 j 6 h 6 y y 6 y y 。 计 算 1 1 n n h h h h 2 2 j1 j j1 j j1 j (8) S(x) M M h 2h 2h h 6 输 出各区间的 s (x

14、)表达式。 S(x)在区间 x ,x ， x ,x （ j 2, ,n-1）上 判 断 x 所在区间 x , x ，并计算插 值 y ， 即S (x 0) S (x 0) 在 x 左右 导数相 等，所以： 3 j1 j j1 j1 在 Matlab 环 境 下根 据 上述 算 法 步骤 进 行编 S (x 0) M M 6 h 3 h y y h function =spline3(X,Y,dY,x0,m) 同样 N=size( X,2); j s0=dY(1); sN=dY(2); h y y h h y y h interval=0.025; j j1 j1 j 3 h 6 6 h 3 h

15、=ze ros(1,N-1) ; j1 j M 2M M for i=1:N-1 h(1,i)=X(i+1) -X(i); e nd h h h h d(1,1)=6*( (Y(1,2) -Y(1,1) )/h( 1,1)- s0)/h(1,1); 1 j j1 j1 j j1 1 d 6 令 j h h /h(1,i-1) )/ ( h( 1,i)+h( 1,i-1) ; end j1 j 兵工 自动 化 自 动测量 与控 制 O. I. Automation 2006 年第 25 卷 第 11 期 Automatic Measurement a nd Contr ol 2006, Vol.

16、 25, No. 11 md( 1,N-1)=1; mu(1,1) =1; M for i=1:N-2 -5.0000 4.0000 4.0000 16.0000 3 2 在区 间 -1.5,0.0 u=h( 1,i+1)/(h( 1,i)+h( 1,i+1) ; mu(1,i+1)=u; md( 1,i)=1- u; end 内，三次 样条函数 S(1) t 2t 2 1 ；在区间 0.0,1.0 p(1,1) =2; q(1,1)=mu( 1,1)/2; for i=2:N-1 p(1,i)=2-md(1,i-1)*q(1,i-1); q(1,i)=mu( 1,i)/p(1,i) ; en

17、d 内，三 次样条函数 S(2) 2 t 三次样条函数 S(3) 2t 4t 2 1 ；在区间 1.0,2.0内， 6t 3 。 p(1,N) =2-md( 1,N-1)*q(1,N- 1) ; y=zeros(1,N); y(1,1)=d(1) /2; x0 所在区 间为 1.0,2.0，函数 在插值点 x0 1.5 for i=2:N y( 1,i)=( d(i)-md( 1,i-1)*y(1,i- 1)/p(1,i) ; end x=z eros(1,N) ; x(1,N) =y(1,N) ; for i=N-1:-1:1 x( 1,i)=y( 1,i)-q( 1,i)*x(1,i+1)

18、 ; end fpr intf (M 为 三 对角 方程 的 解 n); M=x; fpr intf (n); syms t; digits (m) ; for i=1:N-1 pp(i)=M(i)*( X( i+1)-t)3/( 6*h(i) )+M(i+1)*(t-X(i) )3 /( 6*h(i) +(Y(i) -M( i) *h( i) 2/6)*(X(i+1)- t)/h( i)+ ( Y(i+1)-M(i+1)*h(i)2/6) *( t-X(i)/h(i); pp(i)=simplify( pp(i) ); c oe ff=sym2poly(pp(i); if length( c

19、oeff )=4 tt=c oeff (1:3) ; coef f(1:4)=0; coeff (2:4) =tt; end if x0X( i) &x0X(i+1) L=i; y0=coeff (1)*x03+coe ff(2) *x02+coeff (3)*x0+coef f(4); e nd val=X( i):inter val:X( i+1); for k=1:le ngth(val) f val(k)=c oeff (1)*val(k)3+coeff (2)*val(k)2 +coeff (3)*val(k)+coeff (4); e nd if mod(i,2) =1 plot(

20、va l,fval,r+) else plot(val,fval,b.) end hold on clea r va l fva l ans=sym(c oeff ,d); ans=poly2sym( ans,t); fprintf(在 区 间 %f ,%f内 n,X(i),X(i+1) ; fprintf(三 次 样条 函数 S( %d)=,i) ; pretty(a ns); end fpr intf (x0 所在 区 间为 %f,%fn,X(L),X(L+1) ; fpr intf (函 数在 插 值点 x0=%f 的值 为 n,x0) ; y0 程序中： X,Y 为输入 结点， dY

21、为两端点一阶导 数矩阵， x0 为待求插值点， m 为有效数字位数。应 用实例 ：已 知函数 y f(x)的函数值 如表 1。 表 1 y f (x)的 函数 值 x -1.5 0 1 2 y 0. 125 -1 1 9 在区间 -1.5, 2上求三 次样条 函数 S(x)，使它 满足边 界条件 S(-1.5) 0.75， S(2) 14。 解 ：输入有效数字 位数 m 5，待求插值点 x0 1.5。在 Matlab 命令窗口中输入 X -1.5 0 1 2； Y 0.125 -1 1 9； dY 0.75 14； x0 1.5； m 5； spline3(X, Y, dY, x0, m)回车

22、，即可得结果。 78 的值为 y0 3.7500，各区间 段的函数曲线如图 1。 10 8 6 4 2 0 -2 -1 0 1 2 图 1 各区 间 段 的函 数 曲线 5 结论 Mat lab 环境下编写求解三次样 条插值的通用程 序，可直接显示 各区间段三次样条函 数的具体表达 式，计算出已给 点的插值，最后显示 各区间分段曲 线图，为三次样条 插值函数的应用提供简 便方法。 参考文献： 1 易 大 义 , 沈 云 宝 , 李 有 法 . 计 算 方法 M. 浙江 : 浙 江 大 学出 版 社 . 2002, (6): 55-63. 2 谢 文 龙 . 三 次 样 条 函 数 的 构造 方

23、 法 J. 江 南 学 院 学 报 , 2000, ( 6) . 3 张 铮 , 杨 文平 , 石博 强 , 等 . MATLAB 程序 设 计与 实 例 应用 M. 北京 : 中国 铁道 出 版社 , 2003. 来 稿 摘 登 摘登编号： 10 06-1576（ 2006） 11-0 3；收稿日期： 2006-05-2 2 基于 实 值反 向选 择算 法 的 机械 设备 故 障检 测 李 海潮 ， 陈 强 （ 1. 江 西理 工大 学 机电 工 程学 院 赣 州 341000； 2. 北 京科 技 大学 信 息工 程学 院 北京 100083） 摘要 ： 选 择 算法 在实 值 编码 的反

24、 向 选择 机理 中 加入 2 种 进化 学 习机 制。一 将检 测 器集 合远 离 自己 空间 ，二是 移动检 测器，并 最优 化 地将 检测 器 分布 在非 己 空间 。将 其 应用 于机 床齿 轮 箱运 行状 态 检测 ，具 有 很高 的故 障 检测 率。 Machine Equipment Fault Detection Based on Real Numb er Reverse Ch oosing Algor ithm Abstract: The choosing algorithm added two kinds of choosing algorithm into the rev

25、er se choosing me chanism of rea l numbe r c oding. At f ir st, the detec tion set was ke pt f ar from the self-spac e; a t sec ond, the detection was moved a nd distributed optimize d in non-self -spa ce. It wa s use d in mac hine tool gea r box to ca rried out state detection, a nd owne d high fault dete ction rate. mu=zer os(1,N-1); md=zer os( 1,N-1); 0 0 3 1