matlab的应用.docx

1、matlab的应用第12章 三次样条众所周知，使用高阶多项式的插值常常产生病态的结果。目前，有多种消除病态的方法。在这些方法中，三次样条是最常用的一种。在MATLAB中，实现基本的三次样条插值的函数有spline，ppval，mkpp和unmkpp。在这些函数中，仅spline在MATLAB参考指南中有说明。下面几节，将展示在M文件函数中实现三次样条的基本特征。12.1 基本特征在三次样条中，要寻找三次多项式，以逼近每对数据点间的曲线。在样条术语中，这些数据点称之为断点。因为，两点只能决定一条直线，而在两点间的曲线可用无限多的三次多项式近似。因此，为使结果具有唯一性。在三次样条中，增加了三次多

2、项式的约束条件。通过限定每个三次多项式的一阶和二阶导数，使其在断点处相等，就可以较好地确定所有内部三次多项式。此外，近似多项式通过这些断点的斜率和曲率是连续的。然而，第一个和最后一个三次多项式在第一个和最后一个断点以外，没有伴随多项式。因此必须通过其它方法确定其余的约束。最常用的方法，也是函数spline所采用的方法，就是采用非扭结(not-a-knot)条件。这个条件强迫第一个和第二个三次多项式的三阶导数相等。对最后一个和倒数第二个三次多项式也做同样地处理。基于上述描述，人们可能猜想到，寻找三次样条多项式需要求解大量的线性方程。实际上，给定N个断点，就要寻找N-1个三次多项式，每个多项式有4

3、个未知系数。这样，所求解的方程组包含有4*(N-1)个未知数。把每个三次多项式列成特殊形式，并且运用各种约束，通过求解N个具有N个未知系数的方程组，就能确定三次多项式。这样，如果有50个断点，就有50个具有50个未知系数的方程组。幸好，用稀疏矩阵，这些方程式能够简明地列出并求解，这就是函数spline所使用的计算未知系数的方法。12.2 分段多项式在最简单的用法中，spline获取数据x和y以及期望值xi，寻找拟合x和y的三次样条内插多项式，然后，计算这些多项式，对每个xi的值，寻找相应的yi。例如：x=0 : 12;y=tan(pi*x/25);xi=linspace(0, 12);yi=s

4、pline(x, y, xi)plot(x, y, o , xi, yi), title( Spline fit )（见图12.1样条拟合）这种方法适合于只需要一组内插值的情况。不过，如果需要从相同数据集里获取另一组内插值，再次计算三次样条系数是没有意义的。在这种情况下，可以调用仅带前两个参量的spline：图12.1 样条拟合pp=spline(x, y)pp = Columns 1 through 7 10.0000 1.0000 12.0000 0 1.0000 2.0000 3.0000 Columns 8 through 14 4.0000 5.0000 6.0000 7.0000

5、8.0000 9.0000 10.0000 Columns 15 through 21 11.0000 12.0000 4.0000 0.0007 0.0007 0.0010 0.0012 Columns 22 through 28 0.0024 0.0019 0.0116 -0.0083 0.1068 -0.1982 1.4948 Columns 29 through 35 1.4948 -0.0001 0.0020 0.0042 0.0072 0.0109 0.0181 Columns 36 through 42 0.0237 0.0586 0.0336 0.3542 -0.2406 4.

6、2439 0.1257 Columns 43 through 49 0.1276 0.1339 0.1454 0.1635 0.1925 0.2344 0.3167 Columns 50 through 56 0.4089 0.7967 0.9102 4.9136 0 0.1263 0.2568 Columns 57 through 63 0.3959 0.5498 0.7265 0.9391 1.2088 1.5757 2.1251 Columns 64 through 65 3.0777 5.2422当采用这种方式调用时，spline返回一个称之为三次样条的pp形式或分段多项式形式的数组。

7、这个数组包含了对于任意一组所期望的内插值和计算三次样条所必须的全部信息。给定pp形式，函数ppval计算该三次样条。例如，yi=ppval(pp, xi);计算先前计算过的同样的yi。类似地，xi2=linspace(10, 12);yi2=ppval(pp, xi2);运用pp形式，在限定的更细区间10，12内，再次计算该三次样条。xi3=10 : 15yi3=ppval(pp, xi3)yi3 = 3.0777 5.2422 15.8945 44.0038 98.5389 188.4689它表明，可在计算三次多项式所覆盖的区间外，计算三次样条。当数据出现在最后一个断点之后或第一个断点之前时

8、，则分别运用最后一个或第一个三次多项式来寻找内插值。上述给定的三次样条pp形式，存储了断点和多项式系数，以及关于三次样条表示的其它信息。因为，所有信息都被存储在单个向量里，所以这种形式在MATLAB中是一种方便的数据结构。当要计算三次样条表示时，必须把pp形式分解成它的各个表示段。在MATLAB中，通过函数unmkpp完成这一过程。运用上述pp形式，该函数给出如下结果：break, coefs, npolys, ncoefs=unmkpp(pp)breaks = Columns 1 through 12 0 1 2 3 4 5 6 7 8 9 10 11 Column 13 12coefs =

9、 0.0007 -0.0001 0.1257 0 0.0007 0.0020 0.1276 0.1263 0.0010 0.0042 0.1339 0.2568 0.0012 0.0072 0.1454 0.3959 0.0024 0.0109 0.1635 0.5498 0.0019 0.0181 0.1925 0.7265 0.0116 0.0237 0.2344 0.9391 -0.0083 0.0586 0.3167 1.2088 0.1068 0.0336 0.4089 1.5757 -0.1982 0.3542 0.7967 2.1251 1.4948 -0.2406 0.9102

10、 3.0777 1.4948 4.2439 4.9136 5.2422npolys = 12ncoefs = 4这里break是断点，coefs是矩阵，它的第i行是第i个三次多项式，npolys是多项式的数目，ncoefs是每个多项式系数的数目。注意，这种形式非常一般，样条多项式不必是三次。这对于样条的积分和微分是很有益的。给定上述分散形式，函数mkpp恢复了pp形式。pp=mkpp(break, coefs)pp = Columns 1 through 7 10.0000 1.0000 12.0000 0 1.0000 2.0000 3.0000 Columns 8 through 14 4

11、.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 Columns 15 through 21 11.0000 12.0000 4.0000 0.0007 0.0007 0.0010 0.0012 Columns 22 through 28 0.0024 0.0019 0.0116 -0.0083 0.1068 -0.1982 1.4948 Columns 29 through 35 1.4948 -0.0001 0.0020 0.0042 0.0072 0.0109 0.0181 Columns 36 through 42 0.0237 0.05

12、86 0.0336 0.3542 -0.2406 4.2439 0.1257 Columns 43 through 49 0.1276 0.1339 0.1454 0.1635 0.1925 0.2344 0.3167 Columns 50 through 56 0.4089 0.7967 0.9102 4.9136 0 0.1263 0.2568 Columns 57 through 63 0.3959 0.5498 0.7265 0.9391 1.2088 1.5757 2.1251 Columns 64 through 65 3.0777 5.2422因为矩阵coefs的大小确定了npo

13、lys和neofs，所以mkpp不需要npolys和ncoefs去重构pp形式。pp形式的数据结构仅在mkpp中给定为pp=10 1 npolys break(:) ncoefs coefs(:)。前两个元素出现在所有的pp形式中，它们作为确认pp形式向量的一种方法。12.3 积分在大多数情况下，需要知道由三次样条所描述的、自变量为x的函数所包含的面积。也就是，如果这个函数记为y=f(x)，我们感兴趣的是计算： 其中，是s(x1)=0式中的x1是第一个样条的断点。因为s(x)由被连接的三次多项式组成，其中第k个三次多项式为：sk(x)=ak(x-xk)3+ bk(x-xk)2+ ck(x-xk

14、)+dk, xk=x= xk+1并且该函数在区间xk=x= xk+1所含的面积为：三次样条下的面积容易求得为： 式中xk=x= xk+1或者 式中xk=xx=(0: .1: 1)*2*pi;y=sin(x); % create rough datapp=spline(x, y); % pp-form fitting rough datappi=spintgrl(pp); % pp-form of integralxi=linspace(0, 2*pi); % finer points for interpolationyi=ppval(pp, xi); % evaluate curveyyi=

15、ppval(ppi, xi); % evaluate integralplot(x, y, o , xi, yi, xi, yyi, - ) % plot results注意，如图12.2所示，它定性地证明了恒等式：图12.2 函数的积分图形12.4 微分如同人们对样条积分感兴趣一样，一个由样条所描述的函数的微分或斜率也是很有用的。给定第k个三次多项式为：sk(x)=ak(x-xk)3+ bk(x-xk)2+ ck(x-xk)+dk, xk=xx=(0 : .1 : 1)*2*pi; % same data as earliery=sin(x);pp=spline(x, y); % pp-fo

16、rm fitting rough datappd=spderiv(pp); % pp-form of derivativexi=linspace(0, 2*pi); % finer points for interpolationyi=ppval(pp, xi); % evaluate curveyyd=ppval(ppd, xi); % evaluate derivativeplot(x, y, o , xi, yi, xi, yyd, - ) % plot results图12.3 函数的微分图形注意，图12.3定性地证明了恒等式：12.5 小结下面的两个表总结了本章所讨论的样条函数。表1

17、2.1三次样条函数yi=spline(x,y,xi)y=f(x)在xi中各点的三次样条插值pp=spline(x,y)返回y=f(x)的分段多项式的表示yi=ppval(pp,xi)计算xi中各点的分段多项式break,coefs,npolys,ncoefs=unmkpp(pp)分解分段多项式的表示pp=mkpp(break,coefs)形成分段多项式的表示表12.2精通MATLAB的样条函数yi=spintgrl(x,y,xi)在xi各点对y=f(x)积分的三次样条的插值ppi=spintgrl(pp)给定y=f(x)的分段多项式的表示，返回y=f(x)积分的分段多项式的表示。yi=spintgrl(pp,xi)给定y=f(x)的分段多项式的表示，计算点xi的值，寻找 y=f(x)积分的分段多项式的描述。yi=spderiv(x,y,xi)在xi各点对y=f(x)微分的三次样条的插值ppi=spderiv(pp)给定y=f(x)的分段多项式的表示，返回y=f(x)微分的分段多项式的表示。yi=spderiv(pp,xi)给定y=f(x)的分段多项式的表示，计算点xi的值，寻找y=f(x)微分的分段多项式的表示。