MATLAB应用 数据处理.doc

第5章 数据处理51 极值最大值和最小值MATLAB提供的求数据序列的最大值和最小值的函数分别为max和min，两个函数的调用格式和操作过程类似。向量的最大值和最小值求一个向量X的最大值的函数有两种调用格式，分别是：(1) y=max(X)：将向量X的最大值max(X)存入y，如果X中包含复数元素，则按模取最大值。例1 求向量x=-43,72,9,16,23,47的最大值x=-43,72,9,16,23,47y=max(x) %求向量x中的最大值(2) y,k=max(X)：将向量X的最大值max(X)存入y，最大值的序号存入k，如果X中包含复数元素，则按模取最大值。例2 求向量x=-43,72,9,16,23,47的最大值及其该元素的位置x=-43,72,9,16,23,47y,k=max(x) %求向量x中的最大值及其该元素的位置求向量X的最小值的函数是min(X)，用法和max(X)完全相同。例3 求向量x=-43,72,9,16,23,47的最小值及其该元素的位置x=-43,72,9,16,23,47z=min(x) %求向量x中的最小值y,m= min (x) %求向量x中的最小值及其该元素的位置矩阵的最大值和最小值求矩阵A的最大值的函数有3种调用格式，分别是：(1) max(A)：给出一个行向量，向量的第i个元素是矩阵A的第i列上的最大值max(a1i,a2i,ami)。例4 求矩阵A=13, -56, 78; 25, 63, -235; 78, 25, 563; 1, 0, -1各列的最大值及其整个矩阵的的最大值A=13, -56, 78; 25, 63, -235; 78, 25, 563; 1, 0, -1max(A), %求A各列的最大值max(max(A) %相当于求矩阵的的最大值(2) Y,U=max(A)：给出行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。例5 求矩阵A=13, -56, 78; 25, 63, -235; 78, 25, 563; 1, 0, -1各列的最大值及其行号A=13, -56, 78; 25, 63, -235; 78, 25, 563; 1, 0, -1Y,U=max(A)(3) max(A,dim)：dim取1或2。dim取1时，该函数和max(A)完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。例6 用max(A,dim)求矩阵A=13, -56, 78; 25, 63, -235; 78, 25, 563; 1, 0, -1每行及每列的最大值A=13, -56, 78; 25, 63, -235; 78, 25, 563; 1, 0, -1max(A,1) %求A每列的最大值max(A,2) %求A每行的最大值求最小值的函数是min，其用法和max完全相同。例7 求矩阵A=13, -56, 78; 25, 63, -235; 78, 25, 563; 1, 0, -1每行及每列的最小值及整个矩阵的最小值A=13, -56, 78; 25, 63, -235; 78, 25, 563; 1, 0, -1min(A,1) %求A每列的最小值min(A,2) %求A每行的最小值min(min(A) %相当于求矩阵的的最小值 两个向量或矩阵对应元素的比较函数max和min还能对两个同型的向量或矩阵进行比较，调用格式为：(1) U=max(A,B)：A,B是两个同型的向量或矩阵，结果U是与A,B同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。例8 求矩阵A=4, 5, 6; 1, 4, 8及B=1, 7, 5; 4, 5, 7所有同一位置上较大元素构成的新矩阵PA=4, 5, 6; 1, 4, 8B=1, 7, 5; 4, 5, 7P=max(A,B)(2) U=max(A,n)：n是一个标量，结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。例8 求矩阵A=4, 5, 6; 1, 4, 8的所有元素与f=4.5比较后较大者构成的新矩阵P1A=4, 5, 6; 1, 4, 8f=4.5P1=max(A,f)min函数的用法和max完全相同。元素排序sort(X)，它给出一个对向量X中的元素按升序排列的新向量。例9 对向量x=-43,72,9,16,23,47进行升序排列x=-43,72,9,16,23,47sort(x)sort函数也可以对矩阵A的各列或各行重新排序，其调用格式为：Y,k=sort(A,dim)其中dim指明对A的列还是行进行排序。若dim=1，则按列排，可省略；若dim=2，则按行排。Y是排序后的矩阵，而k记录Y中的元素在A中位置。例10 对矩阵A=1, -8, 5; 4, 12, 6; 13, 7, -13进行各种排序A=1, -8, 5; 4, 12, 6; 13, 7, -13,sort(A) %对A的每列按升序排列sort(A,2) %对A的每行按升序排列-sort(-A,2) %对A的每行按降序排列Y,k=sort(A) %对A的按列升序排列，并将每个元素所在的行号送到矩阵k标准方差设有N个数据组成数据序列x1,x2,.xN，这些数据的标准方差为 或 其中，在MATLAB中，提供了计算数据序列的标准方差的函数std。对于向量X，std(X)给出一个标准方差。例11 随机地取7个活塞环，测得它们的直径为(以mm计)74.001, 74.005, 74,003,74.001, 74.000, 73.998, 74.002，求它的标准方差X = 74.001, 74.005, 74,003,74.001, 74.000, 73.998, 74.002std(X)对于矩阵A，std(A)给出一个行向量，它的各个元素便是矩阵A各列或各行的标准方差。std函数的一般调用格式为：Y=std(A,flag,dim)其中dim取1或2。当dim=1时，求各列元素的标准方差；当dim=2时，则求各行元素的标准方差。flag取0或1，当flag=0时，按公式1所列公式计算标准方差，当flag=1时，按公式2所列公式计算标准方差。缺省flag=0，dim=1。例12 对矩阵A=4, 5, 6; 1, 4, 8，从不同维方向求其标准方差A=4, 5, 6; 1, 4, 8, Y1=std(A,0,1); %按按公式1各列元素的标准方差Y2=std(A,1,1); %按按公式2各列元素的标准方差Y3=std(A,0,2); %按按公式1各行元素的标准方差Y4=std(A,1,2); %按按公式2各行元素的标准方差5 . 6 多项式运算在MATLAB中的实现7. 5. 1 多项式及其系数向量1多项式的向量表示法n次多项式表达为：，是n+1项之和。在MATLAB中，多项式是用它的系数构成的行向量来表示的，这个行向量完全等价于它代表的多项式，因此常把“多项式系数向量”就称为“多项式”。该向量的分量自左向右，依次表示多项式高次幂到低次幂项的系数，缺少的幂次项，其系数必须用零填补。如，在MATLAB中多项式，被表示成一个向量p = 3 5 0 -7 9。2多项式转换指令可以将多项式系数向量恢复成多项式的数学形式，即把系数行向量变成易读的表达形式，使用的指令及格式为：poly2str(p, t ) 输入的参数p为多项式系数向量； 输入的t是输出多项式的变量，必须用引号界定，且不得缺省； 输出为一个多项式，其系数是p的各分量，变量为界定的字母t。例如，键入：p = 3 5 0 -7 9; %代表多项式3x4+5x3-7x+9pp = poly2str(p, t )Enter后得到下面的形式：3t4 + 5t3 - 7t + 93特殊多项式的创建一般情况下，创建多项式系数向量需要用键盘输入。但若已知多项式的全部根，则可以用poly函数建立起该多项式，其调用格式为：p=poly(s)若s为具有n个元素的向量，则poly(s)建立以s为其根的多项式，且将该多项式的系数赋给向量p。p与s各分量的排序无关。例如，若键入：s = a b c d, p = poly (s)，这里a , b, c, d都是具体数值时，回车则得出p = (x a)(x b)(x c)(x d)展开成多项式的系数向量。例 5 一 3 求一个多项式，使它的根为3, 5，-7, 9, 0。解 键入：s = 3, 5, -7, 9, 0; p = poly ( s )回车得出： p = 1 -10 -32 474 -945 0再键入：poly2str(p, x )回车得出：ans = x 5 - 10 x 4 - 32 x 3 + 474 x 2 - 945 x这个多项式的数学表示为：。改变输入参量s中各分量的排序，输出向量p不变。例6-22 已知 f(x)=3x5+4x3-5x2-7.2x+5(1) 计算f(x)=0 的全部根。(2) 由方程f(x)=0的根构造一个多项式g(x)，并与f(x)进行对比。P=3,0,4,-5,-7.2,5;X=roots(P) %求方程f(x)=0的根G=poly(X) %求多项式g(x)将这个结果乘以3，就与f(x)一致7. 5. 2 多项式运算在MATLAB中多项式间的运算完全转换为系数向量间的运算，不同的运算对多项式系数向量有一定的要求，在使用中需加注意。1多项式加减与矩阵加减一样，维数相等的系数向量方可进行加减运算，这意味着多项式同幂次对齐后才可加减。当它们的次数相同时，可以直接对多项式的系数向量进行加减运算。当它们的次数不同时，应该把次数低的多项式无高次项部分用0系数表示。也就是必须用“零”填补缺少的高幂次项系数，使它们的系数向量维数相等。例5 一 4 已知，求两个多项式的和解：p1 = 3 2 6; p2 = 3 5 0 -7 9;p1 = 0 0 p1; p = p1+ p2回车得出：p = 3 5 3 -5 15键入： poly2str (p, x)回车得出：ans = 3 x4 + 5 x3 + 3 x2 - 5 x + 15例2 计算 a=1, -2, 5, 3; b=0, 0, 6, -1; c=a+b例3 设，求f(x)+g(x)f=3, -5, 2, -7, 5, 6; g=3, 5, -3; g1=0, 0, 0, g;f+g1, f-g12多项式乘法两个多项式相乘，实际上是求它们系数向量的卷积，指令及使用格式为：conv (p1, p2)输出一个多项式系数向量，其维数是p1和p2维数之和减1。例5 一 5 已知，求这两个多项式的乘积。解：p1 = 3 2 0 -7 2; p2 = 3 -1 0 6; pp =conv (p1, p2)回车得出：pp = 9 3 -2 -3 25 -2 -42 12键入： poly2str(pp, x)回车得出：ans = 9 x7 + 3 x6 - 2 x5 - 3 x4 + 25 x3 - 2 x2 - 42 x + 12表明：3 多项式除法多项式相除就是多项式系数向量的解卷积（卷积的逆运算），指令及使用格式为：Q, r=deconv(p1, p2)表示p1除以p2，给出商式Q(x)，余式r(x)。Q,和r仍为多项式系数向量例5 一 6 利用例5 一 5 中的向量，计算pp / p2和pp / p1。解：q, r = deconv (pp, p1)回车得出：q = 3.0000 -1.0000 0 6.0000r = 1.0e-014 * 0 0 0 0.3553 0 0 -0.7105 0.1776对照例5 一 5 中的数据，可验证其结果，余数数量级很小，仅为10-14，可以忽略。4多项式求导多项式的导函数p=polyder(P)：求多项式P的导函数p=polyder(P,Q)：求PQ的导函数p,q=polyder(P,Q)：求P/Q的导函数，导函数的分子存入p，分母存入q。参数P,Q是多项式的向量表示，p,q也是多项式的向量表示。例4 求有理分式的导函数P=3, 5, 0, -8, 1, -5; %有理分式分子Q=10, 5, 0, 0, 6, 0, 0, 7, -1, 0, -100; %有理分式分母p,q=polyder(P,Q)5多项式求值计算用“数”或“矩阵”代替多项式中变量得出的数值结果时，用下述指令：polyval ( p, x0) 输入参数p为多项式系数向量； 输入参数x0为数值时，得出多项式函数在x0处的值； 输入参数x0为矩阵（或向量）时，输出按数组算法得出的矩阵多项式“值”，即与x0同维的数值矩阵（或向量）。例 5 一 7 已知多项式，求当x = -2，(4, 2, 3)和的取值。解：p1 = 1 0 2 -7 -2; polyval (pl, -2)回车得出： ans = 36键入： polyval (p1 , 4 2 3)回车得出：ans = 258 8 76 键入： polyval (p1 , 2 4 ; 3 -2 )回车得出：ans = 8 25876 3652插值法和数据拟合当变量x在其变化范围内取定一个数值时，相关变量y按一定的法则总有一个或几个定数和它对应，就叫y是x的函数，通常用yF (x)的形式来表示。广义地说， F只表示y与 x 之间的一种对应关系，表达对应关系 F 的具体方法一般有下述3种。 解析法：用数学解析式表示x和y的对应关系，如：y =f(x)x3sinx。 列表法：用一系列的x值和对应的y值列成表，如平方表、三角函数表等。 图示法：在 X - Y 平面上将对应的x和y的数值画成一条曲线。工程技术实践中，特别是在实验测试中，往往只能得到两个相关变量的一系列离散值，它们间的函数关系就只得用列表法或图示法表示。但是，有时因为某种原因希望把这种关系转换成解析表达式。另外，有些变量间的关系虽然可以用解析法表示，由于数学解析式过于繁杂不便使用，也希望用一个简单的解析表达式来近似地代替它。像这样利用简单解析式（x）近似地代替列表法、图示法或复杂解析式表示的函数F(x）一类问题，即寻找一个满足（x）F（x）的简单函数（x）的问题，都可称之为函数逼近。下面介绍的插值法和数据拟合法就是函数逼近的两种重要方法。多项式插值设函数yF(x)是以表5 一 1给出的，即给出了一系列的x和与其对应的y值。表 5 一 1 函数函数yF(x)的列表法表示xx0x1x2xnyy0y1y2yn其中xa , b。 这些数据点反映了一个函数的关系y=F(x)，然而并不知道F(x)的解析式。数值插值的任务是构造一个函数y=（x），使得对于xi(i=0,2,n)有(xi)=F(xi)，而且在两个相邻的采样点(xi, xi+1)，（x）能光滑过渡。如果被插值的函数F(x)是光滑的，而且采样点足够密，一般在采样区间内，F(x)与（x）比较接近。要构造一个简单的解析式（x），它应满足下述的插值原则：（x）F（x）, xa , b;（xi）=yi， i = 0, l , 2 ， ， n 。称（x）为F（x）的插值函数，点x0，x1，x2， ，xn为插值节点（也叫样本点）。由插值原则可知，插值函数在样本点上的取值必须等于已知函数在样本点上的对应值。插值函数（x）可以选用不同的函数形式，如三角多项式、有理函数等。但最常选用的是代数多项式Pn (x) = anxn + an-1xn-1 alx + a0，因为多项式具有形式简单、计算方便、存在各阶导数等良好性质，下边介绍的选用（x）为Pn (x)的代数多项式插值法。代数多项式插值的基本原理1构造n次插值多项式假设对于函数yF(x）给出了一组函数值yiF(xi) (i = 0 , 1 , 2 , 3 ，n）的列表法表示，如表5 一 l 所示。表中x0，x1，x2， ，xn a , b，它们是互异的（n + l）个插值节点。y0，y1，y2， ，yn是分别与这些节点对应的函数值，构造一个不超过n次的多项式： (5 一 l )如果它满足插值原则 （xi）=yi ，即 Pn（xi )yi，( i = 0 , 1 , 2 ， ， n ) ，就被称为F(x）的n次插值多项式。2 确定n次多项式系数求这个多项式就是求出它的（n + l）个系数ai ( i = 0 , 1 , 2 ， ， n）。依次取 x = x0，x1，x2， ，xn，相应地取y0，y1，y2， ，yn，并使其满足插值原则 ：Pn（xi )yi，( i = 0 , 1 , 2 ， ， n )，则可得出系数ai满足的（n + l）个方程构成的方程组： (5 一 2)式（ 5 一 2 ）也可写成下述形式：这个方程组的系数行列式是范德蒙（Vandermonde）行列式： (5 一 3)由范德蒙行列式性质可知，只要xi互异，则V (x0，x1，x2， ，xn） 0 。于是根据克莱姆法则可以知道，方程组（ 5 一 2 ）存在唯一解。也就是说，只要（n + l）个样本点的值x0，x1，x2， ，xn互不相等，原则上就能由式（5 一 2）求出（n + l）个系数 ai，从而就可以得到n次插值多项式Pn（x )。然而，当n很大时，通过解方程组（5 一 2）求系数ai的工作是很烦琐的。为此，需要找出了求解方程组（5 一 2）的简捷方法，下面介绍两种。两种常见插值法求解方程组找出插值多项式系数的方法很多，以拉格朗日（Lagrange）法和牛顿（Newton) 法最为常见，原理简单应用广泛。下面从线性插值、二次插值到 n 次插值，分别对它们予以介绍。1线性插值如果已知函数y F ( x ）列表法表示（表 5 一 2 ) ，表 5 一 2 函数y F ( x ）xx0x1yy0y1这时n = 1 ，用代数插值多项式可写成P1( x ) = a0al x，它满足插值原则 ：Pn（xi )yi ( i = 0 , 1 )时可有：解这个方程组可得出系数：，对a1，和a0进行不同的变换，可得出下面两种形式的插值函数。1 ) Lagrange 线性插值把a0，a1代入插值多项式P1( x ) = a0al x中，变换成直线的两点式表达式：其中，和分别称为节点x0和x1的一次插值基函数。插值函数P1( x )是这两个基函数的线性组合，组合系数就是对应节点上的函数值。这种形式的插值函数叫作Lagrange插值多项式。2 ) Newton 线性插值把a0，a1代入插值多项式P1( x ) = a0al x，变换成直线的点斜式表达式： （55）式中的是函数y = F (x）在xi、xj 处的一阶均差，等于自变量增加一个单位时的函数增量。以均差表示的插值多项式函数（5 一 5) ，称为 Newton 插值多项式。2二次插值如果已知函数y F(x)的列表法表示为表5 一 3，表 5 一 3 函数y F(x)xx0x1x2yy0y1y2则插值多项式最高次数n 2，代数插值多项式为P2(x) = a0al x+a2x2。据插值原则：Pn（xi )yi，( i = 0 , 1 , 2 )可知，系数ai应满足方程组（5 一 2) ，但这里 n = 2，只有三个方程。当各节点的值x0，x1，x2互不相等时，范德蒙（Vandermonde）行列式V (x0，x1，x2） 0，方程组存在唯一解。仿照线性插值的推导方法，有下述两种构成二次插值多项式的方法。 1 ) Lagrange 二次插值可以像推导Lagrange一次插值那样，得出二次插值多项式：P2(x)由二次插值的基函数Li(x), (i = 0, 1, 2）线性组合而成，基函数应该满足：各节点的基函数可根据上式算出：如取i = 0 , L0(x)在x = x1，和x = x2处为零，因此可知L0(x)含有（x 一 xl) (x 一 x2）。又因为L0(x)是一个不多于2次的多项式，不妨设L0(x)= A (x 一 xl) (x 一 x2) ，代入x x0时满足的条件L0(x) = 1，可得 A = l(x0 - xl) (x0 - x2)，于是可以得出L0(x)的表达式。依此方法也可得出L1(x)和L2(x)，把它们一并列在式（5 一 7）中： （5 一 7）于是，把式（5 一 7）代入式（7 一 6）就构成了Lagrange二次插值多项式函数P2(x)。2 ) Newton 二次插值仿照构成Newton一次插值多项式的方法，设二次插值多项式的形式为：满足插值原则 ：Pn（xi )yi，( i = 0 , 1 , 2 )时应有p2(x0）y0，由此可得A = p2(x0）y0。再利用p2(x1）y1，得：称为一阶均差。利用p2(x2）y2，得：这个结果是一阶均差的均差，称为二阶均差。二阶均差的一般形式写成： （5 一 8）于是得到 Newton 二次插值多项式： （5 一 9）3. n次插值1 ) Lagrange 式的 n 次插值若知道函数在（n + l）个节点上的值：( x0，y0）, ( x1，y1）, ( x2，y2），( xn，yn），利用这些数据可以构造出代数插值多项式（5 一 1）。为此，仿照构造二次Lagrange插值多项式的方法，求出相应的n次插值多项式函数： （5 一 10）按照式（5 一 7) ，在节点xi上的Lagrange插值基函数可写成： （5 一 11）若令：，则式（5 一 11）中的，于是n次Lagrange插值多项式（5 一 10）可以写成： （5 一 12）2 ) Newton式的n次插值法仿照前边求出Newton二次插值多项式的办法，首先构造出n阶均差的递推公式： （5 一 13）由式（7 一 13）可以得出在 xi，xj处的一阶均差为：在 xi，xj，xk处的二阶均差为：在 xi，xj，xk，xh处的三阶均差为：以此类推，可以求出F的各阶均差，从而得到 n次牛顿插值多项式： （5 一 14）插值多项式的误差估计用n次插值多项式 近似替代函数F (x）时，设产生的截断误差为Rn(x)，则： （5 一 15）由可知，为 n 次插值多项式的余项。为了估计插值多项式的误差，需要引入余项定理。1余项定理当函数F (x）足够光滑，即满足以下条件。 F (x）在区间a , b上连续。 F (x）具有直至（n + l）阶导数，则对一切xa , b，总存在相应的点a , b，使得： (5 - 16)其中，2 误差估计在式（5 一 16）中，虽然存在，但的值是难以确定的。如果能估算出在区间a，b上的上界 M ，即，则可以得出截断误差的范围： (5 - 17)根据式（5 一 17 ) ，只要找出正数M ，就可以估算出用替代函数 时产生的截断误差 ，实用中往往是根据函数F (x）的具体性质估算M大小的。5 . 3 分段三次插值和三次样条插值仅从截断误差公式来看，用插值多项式近似替代函数时，似乎a , b之间的分点数n越多，插值多项式的次数越高，产生的截断误差就越小，实际上并非如此。 Runge和Faber的工作证明，高次插值多项式并不一定都能收敛到被插值的函数上，而且还增加了许多工作量。例如，将函数 用4次和10次插值多项式函数和替代时，高次插值多项式并不理想，这可从图5 一 1看出。图5 一 1 函数及其插值多项式和的曲线比较图中实线是函数的曲线，两条虚线分别是用插值多项式函数和画的(图中有误，和曲线应该互换)，虽然多项式插值函数都通过了样本点，但是在x0.9时比偏离函数曲线更远，逼近效果更差。5 . 3 . 1 分段三次 Hermite 插值为了利用多项式插值方法而又克服高次插值多项式的缺陷，便引入了分段插值的概念。它的基本思想是把函数在整个区间上分成许多段，每段都选用适当的低次插值多项式代替函数，整体上按一定的光滑性要求连接起来，构成一个分段的插值函数。为此，把函数的自变量x在区间a , b上用（n 一 l）个节点分割成 n 段：a = x0 xl x2 xn-1 xn b 根据这些节点的取值xi、在节点上的函数值和导数值 (i = 0, 1, 2，n )，可以构造一个分段三次插值函数，它满足下述条件。 ， (i = 0, 1, 2，n )。 在每个小区间xi，xi+1（i = 0, 1, 2，n 一 l）上，都是一个三次多项式：把这样构成的分段三次函数称为分段三次Hermite插值函数，它的各小段均为三次多项式，而整体上具有一阶连续导数。5. 3 . 2 三次样条插值的基本原理（“样条”一词源于过去绘图员使用的一种绘图工具样条，它是用富于弹性、能弯曲的木条（或塑料）制成的软尺，把它弯折靠近所有的基点用画笔沿样条就可以画出连接基点的光滑曲线。）三次样条插值也是一种分段插值方法，用分段的三次多项式构造成一个整体上具有函数、一阶和二阶导函数连续的函数，近似地替代已知函数。假设已知函数在区间a , b上的（n + l）个节点a = x0 xl x2 xn-1 xn b及其对应的函数值，( i = 0 , 1 , 2 ， ， n )，即给出（n + l）组样本点数据（x0，y0）, (x1 ，y1)，( xn，yn），可以在a , b上构造一个函数，满足下述条件。 , ( i = 0, 1, 2，n)，即满足插值原则（xi）=yi， i = 0, l , 2 ， ， n 。 在每个小区间xi，xi+1（i = 0, 1, 2，n 一 l）上都是一个三次多项式： ( 5 一 18 ) , 和在a , b上连续。可见，是一个光滑的分段函数，这样的函数称为三次样条（Spline）插值函数。构造的函数是由n个小区间上的分段函数组成，根据条件，每个小区间上构造出一个三次多项式，第 i 个小区间上的三次多项式为，共有n个多项式，每个多项式有4个待定系数。要确定这n个多项式，就需要确定4n个系数 ( i = 0, 1, 2，n - 1）。为此，应该找到包含这些系数的4n个独立方程。根据满足的条件，在所有节点上可得出（n + l）个条件方程：, ( i = 0, 1, 2，n) ( 5 一 19)根据满足的条件，除两端点外在所有节点上，又可得出3(n -1)个方程：, ( i = 0, 1, 2，n) ( 5 一 20)由式（5 一 19）和式（5 一 20）可知共有（4n 一 2）个独立方程，还差两个。通常的办法是在区间a , b的两个端点上各加一个条件，即称之为边界条件。常用的三种边界条件是： 已知和，特别是当取时，称为自然边界条件； 已知和，即已知两端点处切线的斜率； 已知和。这样，在已有的（4n 一 2）个条件方程基础上，再加上任何一种边界条件，即可求出这4n个系数，从而就可以求得三次样条插值函数。5. 3 . 3 三次样条插值函数的一种具体推导方法作为举例，这里介绍一种推导三次样条函数的具体方法。设所求三次样条函数在n个子区间中的任一子区间xi, xi+1（i = 0, 1, 2，n - l）上的三次多项式为：，则可得出：最后一个是线性方程，设，则可以用两点式方程表示：， (5-21)式中，。对式（5 一 21）两边进行一次积分，得出：， (5-22)再对式（5 一 22）两边进行一次积分，得出：， (5-23)式中A , B都是积分常数。式（5 一23 ）中的S (x）应该满足插值原则，即满足式（5 一 19)：，xi，yi是已知的，解这两个方程，得出A （也可先求出B)，代入式(5 一 22）得出：， (5-24)在a , b上连续，所以在相邻两个小区间的分界点xi（节点）上取值相等：，i=1, 2, , n-1 (5-25)由式（5 一 24）和式（5 一 25）便可得出：注：点，点，因此，在代入式（5 一 24）时区间长度分别用和。移项整理可得：若令 b，它们均为常数，于是上式变成： (i=1, 2, , n-1) (5-26)这（n 一 l）个方程中共有（n + l）个未知量：M0，M1，Mn，要求出它们尚缺两个条件方程，这就需要借助边界条件。如果已知M0和Mn，则可由式（5 一 26）得出所有Mi ( i = 0, 1，n) ，把它们和Ai, Bi一并代入式（5 一 23) ，就可求出各子区间上的分段函数。为了再找出两个条件方程，常用一种由内部外推的方法，即在区间a, b起点上用和表示，终点上用和表示。假设，则可得出：又假设，可得到：于是，式（5 一 19）和式（5 一 20）再加上这两个方程，共有4n个条件方程，完全可以确定4n个系数，从而可求出三次插值样条函数。这两个外推条件相当于让两端点处的二阶导数是线性的，即让曲线顺势外延。54 插值法在MATLAB中的实现从已知的一些离散数据点及其函数值，即函数的列表法表示，推求出未知点上函数值的所谓插值方法，在科技工作中应用十分广泛，如查对数表、三角函数表中都会遇到这类插值问题。MATLAB中设有许多插值指令，这里仅介绍最常用的一元函数插值指令，它可以使前面讲过的理论得以用计算机实现。5 . 4 . 1 一元函数的插值（查表）指令interp1该指令的调用格式为：yk = interp1(x, y, xk, method )函数根据数据点x,y的值，计算函数在xk处的值。 输入参数x和y为已知的两个同维向量x =xi1n和 y =yi1n，满足函数yi = F(xi)关系，它们是进行“造表”的根据，把xi, yi称为样本点（插值节点）。 输出量yk是与xk对应的函数值。插值点，可以是数值、向量或矩阵，yk与xk维数相同，其元素一一对应。 用单引号界定的method有 4 种参数可供选择：nearest 最近插值根据已知插值点与已知数据点的远近程度进行插值。插值优先选择较近的数据点进行插值。linear 线性插值把与插值点靠近的两个数据点用直线连接起来，然后在直线上选取对应的插值点的数据。省略method时，即默认为此项。pchip（或cubic）分段三次插值用分段三次多项式Hermite插值（Piecewise Cubic Herlnite Interpolating Polynomial）曲线，依次连接相邻样本点，整体上具有函数及其一阶导数连续性。spline 三次样条插值用分段三次多项式曲线光滑地连接相邻样本点，整体上具有函数、一阶和二阶导数连续性，插值点xk可以在区间外的附近取值，可以是数值、向量或矩阵，yk与xk同维。这个指令并不输出插值多项式函数，只输出插值点上的函数值。这就相当于根据数据对 (x , y) “造表”，然后查出对应于xk的函数值yk，所以又称为查表指令。5. 4. 2 三次插值和三次样条插值指令对于三次插值和三次样条插值，MATLAB中都设有专用指令。三次插值指令调用格式为：yk = pchip (x, y, xk) 三次样条插值指令调用格式为：yk = spline (x, y, xk)参数x, y, xk, yk的意义及要求与interp1中的完全一样，插值效果与interp1中参数method分别选用pchip和spline等价。它允许xk在区间min(x), max(x)外的附近取值，因为程序设计中在两个端点处，都使用了由内部外推的方法。例5 一 1已知y F(x）的函数关系中，当y -0.5 0.3 0 0.8 0.1 0.5 -0.3 0.2 0.3 -0.7 0时，对应的有x =0: length (y) -1。用“最近插值”和“线性插值”法，分别求出当x1 = 0:0.1:length(y)-1时所对应的函数值，并画出曲线。解：题设的x和y是样本点，所以键入：y = -0.5 0.3 0 0.8 0.1 0.5 -0.3 0.2 0.3 -0.7 0; x =0: length (y) -1; %输入样本点x1 = 0:0.1:length (y) -1;yl = interp1 (x, y, x1, nearest); y2 = interp1(x, y, x1, linear); %求出xl的函数值plot(x, y, *, x1, yl, :, x1, y2, r), legend (样本点, 最近插值, 线性插值)如果去掉“求出xl的函数值”前面一行程序末尾的分号，就会输出相应的函数值。直接使用“pchip”和“spline”两个绘图指令，跟用interp1指令时method分别选用pchip和spline是一样的。下面是用三次插值和三次样条作图的例子。例 5 一 2 已知y F(x)函数关系：当y 1 3 0 20 20 4 18时，对应的x =0: length(y）- l。用“三次插值”和“三次样条插值”法，分别求出当 x1 = - 0.5 : 0.2: 5.5时所对应的函数值并画出曲线。解：题中的向量x和y是样本点数据，所以键入：y = 1 3 0 20 20 4 18; x =0: length(y)- 1; x1 = - 0.5 : 0.2: 5.5;y1 = interp1(x, y, x1, pchip);y2 = spline(x, y, x1); plot (x, y, *, x1, y1, r, x1, y2, -.), grid,.legend(样本点,三次插值,三次样条插值);由图可见，用两种方法画出的曲线在样本点之间取值是有差异的，用三次样条插值时曲线的光滑程度要好一些。另外，在限定区间外边界处，三次样条插值和三次插值法的外延趋势不同，差异越发明显。例1 已知的数据点来自，根据生成的数据进行插值处理，生成较平滑的曲线。format long %用15位数字表示结果x=0:.12:1; %生成数据采样点y=(x.2-3*x+5).*exp(-5*x).*sin(x);figure,plot(x,y,x,y,o) %生成数据点x1=0:.02:1; %生成插值采样点y0=(x1.2-3*x1+5).*exp(-5*x1).*sin(x1); %精确值y1=interp1(x,y,x1); %线性插值，缺省项目y2=interp1(x,y,x1,nearest); %最近点插值y3=interp1(x,y,x1,cubic); %3次多项式插值y4=interp1(x,y,x1,spline); %3次样条插值figure,plot(x1,y1, b:, x1,y2, g-., x1,y3, r-, x1,y4, c:,x,y,o,x1,y0) legend(liner,nearest,cubic,spline,数据点,精确值)max(abs(y0(1:49)-y1(1:49), max(abs(y0(1:49)-y2(1:49), max(abs(y0(1:49)-y3(1:49), max(abs(y0(1:49)-y4(1:49) %4种插值的最大误差format short %回到缺省状态很显然，线性插值得到的曲线是很粗糙的；最近点插值得到的曲线更差；采用3次多项式插值和3次样条插值所得到的曲线更接近理论值。例2 概率积分的数据如下表。0.46000.47000.48000.49000.48470.49380.50270.5117用不同的插值方法计算f(0.472)x=0.46:0.01:0.49f=0.4846555,0.4937542,0.5027498,0.5116683format long %用15位数字表示结果interp1(x,f,0.472) interp1(x,f,0.472,nearest)interp1(x,f,0.472,cubic)interp1(x,f,0.472,spline)format short %回到缺省状态5. 5 数据的曲线拟合与数值插值类似，曲线拟合的目的也是用一个比较简单的函数去逼近一个复杂的或未知的函数。数值插值要求逼近函数在采样点与被逼近函数相等。实际上，由于实验或测量误差的存在，所得到的数据本身不一定准确。这时，强求逼近函数通过采样点，显然是不合理的。为此，构造一个函数y=去逼近F(x)，但不要求严格通过采样点。只是希望尽量接近这些采样点。也就是说，使逼近函数与被逼近函数F(x)的在采样点的误差di= - F(xi)在某种意义上最小。寻找“很好逼近”的函数有多种方法，这里只介绍常用的多项式最小二乘法。5. 5. 1 数据曲线拟合的最小二乘法已知两组相关数据x = x1 , x2，xn和y y1 , y2，yn之间存在某种函数关系，如果用一个解析式近似代替，则这个取成m次代数多项式不失为一种最佳选择，因为任何连续函数，至少在一个较小的邻域内都可以用代数多项式任意逼近。若取： ( 5 一 27 )逼近法有许多种途径，下面介绍数据拟合中的最小二乘法。最小二乘法拟合是要求在所有n个样本点xj (j=1, 2, , n)处，多项式取值与函数值yj偏差的平方之和达到最小，也就是使“很好地逼近”函数值，为此，要求下面表达式中的R达到最小（在插值法中要求，即R = 0，且一般样本点数 n = m + 1 )： (5一 28)这种要求更符合实际需要，因为“偏差的平方和”尽可能小就保证了偏差绝对值尽可能小，这正是对实测数据的希望。由式（5 一 28）可知， R是待求变量 a0，a1, a2，am的函数：R = R (a0，a1, a2，am）。使R尽可能地小，就归结为求多元函数R的极小值。这可以用数学分析中求极值的方法，即让R对ai ( i = 0, 1, 2，m）的偏导数都等于零，从而求出各项系数a0，a1, a2，am满足的方程为：移项可得： (5 一 29)这是由m个方程构成的方程组，即下述的矩阵方程：式中的均指j从1到n取和。求解方程组（5 一 29) ，就可得出（m + l）个系数a0，a1, a2，am，也就求得了拟合多项式。下面介绍求解这个方程组的近似解的方法。5 . 5 . 2 超定方程组的最小二