matlab的应用.docx
《matlab的应用.docx》由会员分享，可在线阅读，更多相关《matlab的应用.docx（11页珍藏版）》请在冰豆网上搜索。
matlab的应用
第12章三次样条
众所周知，使用高阶多项式的插值常常产生病态的结果。
目前，有多种消除病态的方法。
在这些方法中，三次样条是最常用的一种。
在MATLAB中，实现基本的三次样条插值的函数有spline，ppval，mkpp和unmkpp。
在这些函数中，仅spline在《MATLAB参考指南》中有说明。
下面几节，将展示在M文件函数中实现三次样条的基本特征。
12.1基本特征
在三次样条中，要寻找三次多项式，以逼近每对数据点间的曲线。
在样条术语中，这些数据点称之为断点。
因为，两点只能决定一条直线，而在两点间的曲线可用无限多的三次多项式近似。
因此，为使结果具有唯一性。
在三次样条中，增加了三次多项式的约束条件。
通过限定每个三次多项式的一阶和二阶导数，使其在断点处相等，就可以较好地确定所有内部三次多项式。
此外，近似多项式通过这些断点的斜率和曲率是连续的。
然而，第一个和最后一个三次多项式在第一个和最后一个断点以外，没有伴随多项式。
因此必须通过其它方法确定其余的约束。
最常用的方法，也是函数spline所采用的方法，就是采用非扭结（not-a-knot）条件。
这个条件强迫第一个和第二个三次多项式的三阶导数相等。
对最后一个和倒数第二个三次多项式也做同样地处理。
基于上述描述，人们可能猜想到，寻找三次样条多项式需要求解大量的线性方程。
实际上，给定N个断点，就要寻找N-1个三次多项式，每个多项式有4个未知系数。
这样，所求解的方程组包含有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=spline（x,y,xi）
>>plot（x,y,‘o‘,xi,yi）,title（‘Splinefit‘）
（见图12.1样条拟合）
这种方法适合于只需要一组内插值的情况。
不过，如果需要从相同数据集里获取另一组内插值，再次计算三次样条系数是没有意义的。
在这种情况下，可以调用仅带前两个参量的spline：
图12.1样条拟合
>>pp=spline（x,y）
pp=
Columns1through7
10.00001.000012.000001.00002.00003.0000
Columns8through14
4.00005.00006.00007.00008.00009.000010.0000
Columns15through21
11.000012.00004.00000.00070.00070.00100.0012
Columns22through28
0.00240.00190.0116-0.00830.1068-0.19821.4948
Columns29through35
1.4948-0.00010.00200.00420.00720.01090.0181
Columns36through42
0.02370.05860.03360.3542-0.24064.24390.1257
Columns43through49
0.12760.13390.14540.16350.19250.23440.3167
Columns50through56
0.40890.79670.91024.913600.12630.2568
Columns57through63
0.39590.54980.72650.93911.20881.57572.1251
Columns64through65
3.07775.2422
当采用这种方式调用时，spline返回一个称之为三次样条的pp形式或分段多项式形式的数组。
这个数组包含了对于任意一组所期望的内插值和计算三次样条所必须的全部信息。
给定pp形式，函数ppval计算该三次样条。
例如，
>>yi=ppval（pp,xi）;
计算先前计算过的同样的yi。
类似地，
>>xi2=linspace（10,12）;
>>yi2=ppval（pp,xi2）;
运用pp形式，在限定的更细区间[10，12]内，再次计算该三次样条。
>>xi3=10:
15
>>yi3=ppval（pp,xi3）
yi3=
3.07775.242215.894544.003898.5389188.4689
它表明，可在计算三次多项式所覆盖的区间外，计算三次样条。
当数据出现在最后一个断点之后或第一个断点之前时，则分别运用最后一个或第一个三次多项式来寻找内插值。
上述给定的三次样条pp形式，存储了断点和多项式系数，以及关于三次样条表示的其它信息。
因为，所有信息都被存储在单个向量里，所以这种形式在MATLAB中是一种方便的数据结构。
当要计算三次样条表示时，必须把pp形式分解成它的各个表示段。
在MATLAB中，通过函数unmkpp完成这一过程。
运用上述pp形式，该函数给出如下结果：
>>[break,coefs,npolys,ncoefs]=unmkpp（pp）
breaks=
Columns1through12
01234567891011
Column13
12
coefs=
0.0007-0.00010.12570
0.00070.00200.12760.1263
0.00100.00420.13390.2568
0.00120.00720.14540.3959
0.00240.01090.16350.5498
0.00190.01810.19250.7265
0.01160.02370.23440.9391
-0.00830.05860.31671.2088
0.10680.03360.40891.5757
-0.19820.35420.79672.1251
1.4948-0.24060.91023.0777
1.49484.24394.91365.2422
npolys=
12
ncoefs=
4
这里break是断点，coefs是矩阵，它的第i行是第i个三次多项式，npolys是多项式的数目，ncoefs是每个多项式系数的数目。
注意，这种形式非常一般，样条多项式不必是三次。
这对于样条的积分和微分是很有益的。
给定上述分散形式，函数mkpp恢复了pp形式。
>>pp=mkpp（break,coefs）
pp=
Columns1through7
10.00001.000012.000001.00002.00003.0000
Columns8through14
4.00005.00006.00007.00008.00009.000010.0000
Columns15through21
11.000012.00004.00000.00070.00070.00100.0012
Columns22through28
0.00240.00190.0116-0.00830.1068-0.19821.4948
Columns29through35
1.4948-0.00010.00200.00420.00720.01090.0181
Columns36through42
0.02370.05860.03360.3542-0.24064.24390.1257
Columns43through49
0.12760.13390.14540.16350.19250.23440.3167
Columns50through56
0.40890.79670.91024.913600.12630.2568
Columns57through63
0.39590.54980.72650.93911.20881.57572.1251
Columns64through65
3.07775.2422
因为矩阵coefs的大小确定了npolys和neofs，所以mkpp不需要npolys和ncoefs去重构pp形式。
pp形式的数据结构仅在mkpp中给定为pp=[101npolysbreak（:
）‘ncoefscoefs（:
）‘]。
前两个元素出现在所有的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）+dk,xk<=x<=xk+1
并且该函数在区间xk<=x<=xk+1所含的面积为：
三次样条下的面积容易求得为：
式中xk<=x<=xk+1
或者
式中xk<=x<=xk+1
上式相加项是所有处理的三次多项式面积的累加和。
照此，因为Sk（x）是一个多项式，所以该面积很容易计算，且在描述S（x）的多项式中可形成常数项。
有了上述理解，积分本身可写成样条形式。
在这种情况下，因为单个多项式具有4阶，所以积分是四次样条。
MATLAB中使用的pp形式支持任意阶的样条，所以在《精通MATLAB工具箱》中的函数spintgrl体现了上述样条积分。
该函数的主体如下：
functionz=spintgrl（x,y,xi）
%SPINTGRLCubicSplineIntegralInterpolation
%YI=SPINTRGL（X,Y,XI）usescubicsplineinterpolationtofitthedateinXandY,integrates
%thesplineandreturnsvaluesoftheintegralevaluatedatthepointsinXI.
%
%PPI=SPINTGRL（PP）returnsthepiecewisepolynomialvectorPPI
%describingtheintegralofthecubicsplinedescribedby
%thepiecewisepolynomialinPP.PPisreturnedbythefunction
%SPLINEandisadatavectorcontainingallinformationto
%evaluateandmanipulateaspline.
%
%YI=SPINTGRL（PP,XI）integratesthecubicsplinesgivenby
%thepiecewisepolynomialPP,andreturnsthevaluesofthe
%integralevaluatedatthepointsinXI.
%
%SeealsoSPLINE,PPVAL,MKPP,UNMKPP,SPDERIV
%Copyright（c）1996byPrentice-Hall,Inc.
ifnargin==3
pp=spline（x,y）
else
pp=x;
end
[br,co,npy,nco]=unmkpp（pp）;%takeapartpp
ifpp
（1）==10
error（‘Splinedatadoesnothavethecorrectform.‘）
end
sf=nco:
-1:
1;%scalefactorsforintegration
ico=[co./sf（ones（npy,1）,:
）zeros（npy,1）];%integralcoefficients
nco=nco+1;%splineorderincreases
fork=2:
npy%findconstantterms
ico（k,nco）=polyval（ico（k-1,:
）,br（k）-br（k-1））;
end
ppi=mkpp（br,ico）;%buildppformforintegral
ifnargin==1
z=ppi;
elseifnargin==2
z=ppval（ppi,y）;
else
z=ppval（ppi,xi）;
end
考虑如下运用spintgrl的例子：
>>x=（0:
.1:
1）*2*pi;
>>y=sin（x）;%createroughdata
>>pp=spline（x,y）;%pp-formfittingroughdata
>>ppi=spintgrl（pp）;%pp-formofintegral
>>xi=linspace（0,2*pi）;%finerpointsforinterpolation
>>yi=ppval（pp,xi）;%evaluatecurve
>>yyi=ppval（ppi,xi）;%evaluateintegral
>>plot（x,y,‘o‘,xi,yi,xi,yyi,‘-‘）%plotresults
注意，如图12.2所示，它定性地证明了恒等式：
图12.2函数的积分图形
12.4微分
如同人们对样条积分感兴趣一样，一个由样条所描述的函数的微分或斜率也是很有用的。
给定第k个三次多项式为：
sk（x）=ak（x-xk）3+bk（x-xk）2+ck（x-xk）+dk,xk<=x<=xk+1
容易写出其导数为：
如同样条积分，样条微分也是一种样条。
不过，在这种情况下，这个多项式为二阶，所以它是二次样条。
基于上述描述，《精通MATLAB工具箱》中的函数spderiv实施样条微分。
spderiv的主体为：
functionz=spderiv（x,y,xi）
%SPDERIVCubicSplineDerivativeInterpolation
%YI=SPDERIV（X,Y,XI）usescubicsplineinterpolationtofitthe
%datainXandY,differentiatesthesplineandreturnsvalues
%ofthesplinederivativesevaluatedatthepointsinXI.
%
%PPD=SPDERIV（PP）returnsthepiecewisepolynomialvectorPPD
%describingthecubicsplinederivativeofthecurvedescribedby
%thepiecewisepolynomialinPP.PPisreturnedbythefunction
%SPLINEandisadatavectorcontainingallinformationto
%evaluateandmanipulateaspline.
%
%YI=SPDERIV（PP,XI）differentiatesthecubicsplinegivenby
%thepiecewisepolynomialPP,andreturnsthevalueofthe
%splinederivativesevaluatedatthepointsinXI.
%
%SeealsoSPLINE,PPVAL,MKPP,UNMKPP,SPINTGRL
%Copyright（c）1996byPrentice-Hall,Inc.
ifnargin==3
pp=spline（x,y）;
else
pp=x;
end
[br,co,npy,nco]=unmkpp（pp）;%takeapartpp
ifnco==1|pp
（1）~=10
error（‘SplinedatadoesnothavethecorrectPPform.‘）
end
sf=nco-1:
-1:
1;%scalefactorsfordifferentiation
dco=sf（ones（npy,1）,:
）.*co（:
1:
nco-1）;%derivativecoefficients
ppd=mkpp（br,dco）;%buildppformforderivative
ifnargin==1
z=ppd;
elseifnargin==2
z=ppval（ppd,y）;
else
z=ppval（ppd,xi）;
end
为演示spderiv的使用，考虑如下例子：
>>x=（0:
.1:
1）*2*pi;%samedataasearlier
>>y=sin（x）;
>>pp=spline（x,y）;%pp-formfittingroughdata
>>ppd=spderiv（pp）;%pp-formofderivative
>>xi=linspace（0,2*pi）;%finerpointsforinterpolation
>>yi=ppval（pp,xi）;%evaluatecurve
>>yyd=ppval（ppd,xi）;%evaluatederivative
>>plot（x,y,‘o‘,xi,yi,xi,yyd,‘-‘）%plotresults
图12.3函数的微分图形
注意，图12.3定性地证明了恒等式：
12.5小结
下面的两个表总结了本章所讨论的样条函数。
表12.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）微分的分段多项式的表示。
升级会员 