三次样条插值算法详解.ppt

* 从而推导出了三次样条S(x)在第k个小区间[xk,xk+1]上的表达式为: 它的系数都是用二阶导数与函数值表示! * 对所有中间节点xk,k=1,2,…,n-1,左边小区间与右边小区间上的三次多项式的一阶导数应当连续! 确定二阶导数 * 三弯矩法基本方程 注意到这个基本方程只包括了n-1个方程!但却有n个二阶导数需要待定,这是一个欠定方程组,还需要根据边界条件再确定两个方程! * 曲率调整样条 这种样条的边界条件是已知两端点的二阶导数值! 这样从三弯矩基本方程可以导数确定其它n-2个待定参数的方程组: * 自然样条 这种样条的边界条件是:已知两端点的二阶导数值为0! 这样从三弯矩基本方程可以导数确定其它n-2个待定参数的方程组: * 固支样条 这种样条的边界条件是:已知两端点的一阶导数值! 根据前面推导过程中得到的样条函数S(x)的一阶导数的表达式(2.11),得方程 * 固支样条 这样从三弯矩基本方程可以导数确定n个待定参数的方程组: * 非扭结样条 这种样条的边界条件是:要求样条S(x)在开始的两个小区间[x0,x1],[x1,x2]上的三阶导数相同,在最后两个小区间[xn-2,xn-1],[xn-1,xn]上的三阶导数相同. 对表达式(2.9)再求一次导数得方程 * 非扭结样条 再由三弯矩基本方程,可得 * 周期样条 这种样条的边界条件是:要求样条S(x)及其导数是以区间长度xn-x0为周期的函数即 这些条件可以确定如下两个方程: * 再由三弯矩基本方程,可得 周期样条 * 在Matlab中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。这个条件强迫第1个和第2个三次多项式的三阶导数相等，对最后一个和倒数第2个三次多项式也做同样地处理。 Matlab中三次样条插值也有现成的函数： y=interp1(x0,y0,x,spline)； y=spline(x0,y0,x)； pp=csape(x0,y0,conds), pp=csape(x0,y0,conds,valconds)， y=ppval(pp,x)。 其中x0,y0是已知数据点，x是插值点，y是插值点的函数值。 对于三次样条插值，我们提倡使用函数csape，csape的返回值是pp形式，要求插值点的近似函数值，必须调用函数ppval。 MATLAB中三次样条函数法实现 * pp=csape(x0,y0,conds,valconds) conds指定插值的边界条件，其值可为： complete 边界为一阶导数，一阶导数的值在valconds参数中给出。 not-a-knot 非扭结条件 periodic 周期条件 second 边界为二阶导数，二阶导数的值在valconds参数中给出，若忽略valconds参数，二阶导数的缺省值为[0, 0]。 MATLAB中三次样条函数法实现 * 例2：第一边界条件的例题 x 1 2 4 5 y 1 3 4 2 y’ 17/8 -19/8 clear; x=[1,2,4,5]; y=[1,3,4,2]; pp=csape(x,y,complete,[17/8,-19/8]); pp.coefs MATLAB代码 * 第一边界条件的例题 pp = form: pp breaks: [1 2 4 5] coefs: [3x4 double] pieces: 3 order: 4 dim: 1 pp.coefs -0.1250 0 2.1250 1.0000 -0.1250 -0.3750 1.7500 3.0000 0.3750 -1.1250 -1.2500 4.0000 分段多项式结构与系数矩阵 * 如下代码求解上述样条问题： x -4 -3 -2 -1 0 1 2 3 4 y 0 0.15 1.12 2.36 2.36 1.46 0.49 0.06 0 y’ 0 0 例3：自然边界条件例题 * x=[-4 -3 -2 -1 0 1 2 3 4]; y=[0 0.15 1.12 2.36 2.36 1.46 0.49 0.06 0]; pp=csape(x,y,second); xx=-4:0.01:4; yy=ppval(pp,xx); hold on; plot(x,y,ok); plot(xx,yy,k-); hold off; MATLAB程序 * pp = form: pp breaks: [-4 -3 -2 -1 0 1 2 3 4] coe