02 数值计算与数据分析.pdf

MATLAB6.0 数学手册 62第第 2 章章 数值计算与数据分析数值计算与数据分析 2.1 基本数学函数基本数学函数 2.1.1 三角函数与双曲函数三角函数与双曲函数 函数 sin、sinh 功能 正弦函数与双曲正弦函数 格式 Y=sin(X)%计算参量 X（可以是向量、矩阵，元素可以是复数）中每一个角度分量的正弦值 Y，所有分量的角度单位为弧度。Y=sinh(X)%计算参量 X 的双曲正弦值 Y 注意：sin(pi)并不是零，而是与浮点精度有关的无穷小量 eps，因为 pi 仅仅是精确值浮点近似的表示值而已；对于复数 Z=x+iy，函数的定义为：sin(x+iy)=sin(x)*cos(y)+i*cos(x)*sin(y)，2ee)zsin(iziz=，2ee)zsin(zz=例 2-1 x=-pi:0.01:pi;plot(x,sin(x)x=-5:0.01:5;plot(x,sinh(x)图形结果为图 2-1。图 2-1 正弦函数与双曲正弦函数图 函数 asin、asinh 功能 反正弦函数与反双曲正弦函数 格式 Y=asin(X)%返回参量 X（可以是向量、矩阵）中每一个元素的反正弦函数值Y。若 X 中有的分量处于-1,1之间，则 Y=asin(X)对应的分量处于-/2,/2之间，若 X 中有分量在区间-1,1之外，则 Y=asin(X)对应的分量为复数。Y=asinh(X)%返回参量 X 中每一个元素的反双曲正弦函数值 Y 第 2 章 数值计算与数据分析 63说明 反正弦函数与反双曲正弦函数的定义为：)z1ziln(izsina2+=，)z1zln(zsinha2+=例 2-2 x=-1:.01:1;plot(x,asin(x)x=-5:.01:5;plot(x,asinh(x)图形结果为图 2-2。图 2-2 反正弦函数与反双曲正弦函数图 函数 cos、cosh 功能 余弦函数与双曲余弦函数 格式 Y=cos(X)%计算参量 X（可以是向量、矩阵，元素可以是复数）中每一个角度分量的余弦值 Y，所有角度分量的单位为弧度。我们要指出的是，cos(pi/2)并不是精确的零，而是与浮点精度有关的无穷小量 eps，因为 pi 仅仅是精确值浮点近似的表示值而已。Y=sinh(X)%计算参量 X 的双曲余弦值 Y 说明 若 X 为复数 z=x+iy，则函数定义为：cos(x+iy)=cos(x)*cos(y)+i*sin(x)*sin(y)，2eezcosiziz+=，2eezcoshzz+=例 2-3 x=-pi:0.01:pi;plot(x,cos(x)x=-5:0.01:5;plot(x,cosh(x)图形结果为图 2-3。图 2-3 余弦函数与双曲余弦函数图 函数 acos、acosh 功能 反余弦函数与反双曲余弦函数 MATLAB6.0 数学手册 64格式 Y=acos(X)%返回参量 X（可以是向量、矩阵）中每一个元素的反余弦函数值 Y。若 X 中有的分量处于-1,1之间，则 Y=acos(X)对应的分量处于0,之间，若X中有分量在区间-1,1之外，则Y=acos(X)对应的分量为复数。Y=asinh(X)%返回参量 X 中每一个元素的反双曲余弦函数 Y 说明 反余弦函数与反双曲余弦函数定义为：)z1iziln(izcosa2+=，)1zzln(zcosha2+=例 2-4 x=-1:.01:1;plot(x,acos(x)x=-5:.01:5;plot(x,acosh(x)图形结果为图 2-4。图 2-4 反余弦函数与反双曲余弦函数图 函数 tan、tanh 功能 正切函数与双曲正切函数 格式 Y=tan(X)%计算参量 X（可以是向量、矩阵，元素可以是复数）中每一个角度分量的正切值 Y，所有角度分量的单位为弧度。我们要指出的是，tan(pi/2)并不是精确的零，而是与浮点精度有关的无穷小量 eps，因为 pi 仅仅是精确值浮点近似的表示值而已。Y=tanh(X)%返回参量 X 中每一个元素的双曲正切函数值 Y 例 2-5 x=(-pi/2)+0.01:0.01:(pi/2)-0.01;%稍微缩小定义域 plot(x,tan(x)x=-5:0.01:5;plot(x,tanh(x)图形结果为图 2-5。图 2-5 正切函数与双曲正切函数图 第 2 章 数值计算与数据分析 65函数 atan、atanh 功能 反正切函数与反双曲正切函数 格式 Y=atan(X)%返回参量 X（可以是向量、矩阵）中每一个元素的反正切函数值 Y。若 X 中有的分量为实数，则 Y=atan(X)对应的分量处于-/2,/2之间。Y=atanh(X)%返回参量 X 中每一个元素的反双曲正切函数值 Y。说明 反正切函数与反双曲正切函数定义为：ziziln2iztana+=，z1z1ln21ztanha+=例 2-6 x=-20:0.01:20;plot(x,atan(x)x=-0.99:0.01:0.99;plot(x,atanh(x)图形结果为图 2-6。图 2-6 反正切函数与反双曲正切函数图 函数 cot、coth 功能 余切函数与双曲余切函数 格式 Y=cot(X)%计算参量 X（可以是向量、矩阵，元素可以是复数）中每一个角度分量的余切值 Y，所有角度分量的单位为弧度。Y=coth(X)%返回参量 X 中每一个元素的双曲余切函数值 Y 例 2-7 x1=-pi+0.01:0.01:-0.01;%去掉奇点 x=0 x2=0.01:0.01:pi-0.01;%做法同上 plot(x1,cot(x1),x2,cot(x2)plot(x1,coth(x1),x2,coth(x2)图形结果为图 2-7。图 2-7 余切函数与双曲余切函数图 MATLAB6.0 数学手册 66函数 acot、acoth 功能 反余切函数与反双曲余切函数 格式 Y=acot(X)%返回参量 X（可以是向量、矩阵）中每一个元素的反余切函数 Y Y=acoth(X)%返回参量 X 中每一个元素的反双曲余切函数值 Y 例 2-8 x1=-2*pi:pi/30:-0.1;x2=0.1:pi/30:2*pi;%去掉奇异点 x=0 plot(x1,acot(x1),x2,acot(x2)x1=-30:0.1:-1.1;x2=1.1:0.1:30;plot(x1,acoth(x1),x2,acoth(x2)图形结果为图 2-8。图 2-8 反余切函数与反双曲余切函数图 函数 sec、sech 功能 正割函数与双曲正割函数 格式 Y=sec(X)%计算参量 X（可以是向量、矩阵，元素可以是复数）中每一个角度分量的正割函数值 Y，所有角度分量的单位为弧度。我们要指出的是，sec(pi/2)并不是无穷大，而是与浮点精度有关的无穷小量 eps的倒数，因为 pi 仅仅是精确值浮点近似的表示值而已。Y=sech(X)%返回参量 X 中每一个元素的双曲正割函数值 Y 例 2-9 x1=-pi/2+0.01:0.01:pi/2-0.01;%去掉奇异点 x=pi/2 x2=pi/2+0.01:0.01:(3*pi/2)-0.01;plot(x1,sec(x1),x2,sec(x2)x=-2*pi:0.01:2*pi;plot(x,sech(x)图形结果为图 2-9。图 2-9 正割函数与双曲正割函数图 第 2 章 数值计算与数据分析 67函数 asec、asech 功能 反正割函数与反双曲正割函数 格式 Y=asec(X)%返回参量 X（可以是向量、矩阵）中每一个元素的反正割函数值Y Y=asech(X)%返回参量 X 中每一个元素的反双曲正割函数值 Y 例 2-10 x1=-5:0.01:-1;x2=1:0.01:5;plot(x1,asec(x1),x2,asec(x2)x=0.01:0.001:1;plot(x,asech(x)图形结果为图 2-10。图 2-10 反正割函数与反双曲正割函数图 函数 csc、csch 功能 余割函数与双曲余割函数 格式 Y=csc(X)%计算参量 X（可以是向量、矩阵，元素可以是复数）中每一个角度分量的余割函数值 Y，所有角度分量的单位为弧度。Y=csch(X)%返回参量 X 中每一个元素的双曲余割函数值 Y 例 2-11 x1=-pi+0.01:0.01:-0.01;x2=0.01:0.01:pi-0.01;%去掉奇异点 x=0 plot(x1,csc(x1),x2,csc(x2)plot(x1,csch(x1),x2,csch(x2)图形结果为图 2-11。图 2-11 余割函数与双曲余割函数图 函数 acsc、acsch MATLAB6.0 数学手册 68功能 反余割函数与反双曲余割函数。格式 Y=asec(X)%返回参量 X（可以是向量、矩阵）中每一个元素的反余割函数值 Y Y=asech(X)%返回参量 X 中每一个元素的反双曲余割函数值 Y 例 2-12 x1=-10:0.01:-1.01;x2=1.01:0.01:10;%去掉奇异点 x=1 plot(x1,acsc(x1),x2,acsc(x2)x1=-20:0.01:-1;x2=1:0.01:20;plot(x1,acsch(x1),x2,acsch(x2)图形结果为图 2-12。图 2-12 反余割函数与反双曲余割函数图 函数 atan2 功能 四象限的反正切函数 格式 P=atan2(Y,X)%返回一与参量 X 和 Y 同型的、与 X 和 Y 元素的实数部分对应的、元素对元素的四象限的反正切函数阵列 P，其中 X 和 Y的虚数部分将忽略。阵列 P 中的元素分布在闭区间-pi,pi上。特定的象限将取决于 sign(Y)与 sign(X)。例 2-13 z=1+2i;r=abs(z);theta=atan2(imag(z),real(z)z=r*exp(i*theta)feather(z);hold on t=0:0.1:2*pi;x=1+sqrt(5)*cos(t);y=sqrt(5)*sin(t);plot(x,y);axis equal;hold off 计算结果为：theta=1.1071 z=1.0000+2.0000i 图形结果为图 2-13。xy-x0y0 x0y0 x0 x0yA=-1.9,-0.2,3.1415926,5.6,7.0,2.4+3.6i;B=fix(A)计算结果为：B=Columns 1 through 4 -1.0000 0 3.0000 5.0000 Columns 5 through 6 7.0000 2.0000+3.0000i 函数 roud 功能 朝最近的方向取整。格式 Y=round(X)%对 X 的每一个元素朝最近的方向取整数部分，返回与 X 同维的数组。对于复数参量 X，则返回一复数，其分量的实数与虚数部分分别取原复数的、朝最近方向的整数部分。例 2-15 A=-1.9,-0.2,3.1415926,5.6,7.0,2.4+3.6i;Y=round(A)计算结果为：Y=Columns 1 through 4 -2.0000 0 3.0000 6.0000 Columns 5 through 6 7.0000 2.0000+4.0000i 函数 floor 功能 朝负无穷大方向取整 格式 B=floor(A)%对 A 的每一个元素朝负无穷大的方向取整数部分，返回与 A 同维的数组。对于复数参量 A，则返回一复数，其分量的实数与虚MATLAB6.0 数学手册 70数部分分别取原复数的、朝负无穷大方向的整数部分。例 2-16 A=-1.9,-0.2,3.1415926,5.6,7.0,2.4+3.6i;F=floor(A)计算结果为：F=Columns 1 through 4 -2.0000 -1.0000 3.0000 5.0000 Columns 5 through 6 7.0000 2.0000+3.0000i 函数 rem 功能 求作除法后的剩余数 格式 R=rem(X,Y)%返回结果 X-fix(X./Y).*Y，其中 X、Y 应为正数。若 X、Y 为浮点数，由于计算机对浮点数的表示的不精确性，则结果将可能是不可意料的。fix(X./Y)为商数 X./Y 朝零方向取的整数部分。若 X与 Y为同符号的，则rem(X,Y)返回的结果与 mod(X,Y)相同，不然，若 X 为正数，则 rem(-X,Y)=mod(-X,Y)-Y。该命令返回的结果在区间0，sign(X)*abs(Y)，若 Y 中有零分量，则相应地返回 NaN。例 2-17 X=12 23 34 45;Y=3 7 2 6;R=rem(X,Y)计算结果为：R=0 2 0 3 函数 ceil 功能 朝正无穷大方向取整 格式 B=floor(A)%对 A 的每一个元素朝正无穷大的方向取整数部分，返回与 A同维的数组。对于复数参量 A，则返回一复数，其分量的实数与虚数部分分别取原复数的、朝正无穷大方向的整数部分。例 2-18 A=-1.9,-0.2,3.1415926,5.6,7.0,2.4+3.6i;B=ceil(A)计算结果为：B=Columns 1 through 4 -1.0000 0 4.0000 6.0000 Columns 5 through 6 7.0000 3.0000+4.0000i 函数 exp 功能 以 e 为底数的指数函数 格式 Y=exp(X)%对参量 X 的每一分量，求以 e 为底数的指数函数 Y。X 中的分量可以为复数。对于复数分量如，z=x+i*y，则相应地计算：ez=ex*(cos(y)+i*sin(y)。第 2 章 数值计算与数据分析 71例 2-19 A=-1.9,-0.2,3.1415926,5.6,7.0,2.4+3.6i;Y=exp(A)计算结果为：Y=1.0e+003*Columns 1 through 4 0.0001 0.0008 0.0231 0.2704 Columns 5 through 6 1.0966 -0.0099-0.0049i 函数 expm 功能 求矩阵的以 e 为底数的指数函数 格式 Y=expm(X)%计算以 e 为底数、x 的每一个元素为指数的指数函数值。若矩阵 x 有小于等于零的特征值，则返回复数的结果。说明 该函数为一内建函数，它有三种计算算法：（1）使用文件 expm1.m 中的用比例法与二次幂算法得到的 Pad 近似值；（2）使用 Taylor 级数近似展开式计算，这种计算在文件 expm2.m 中。但这种一般计算方法是不可取的，通常计算是缓慢且不精确的；（3）在文件 expm3.m 中，先是将矩阵对角线化，再把函数计算出相应的的特征向量，最后转换过来。但当输入的矩阵没有与矩阵阶数相同的特征向量个数时，就会出现错误。例 2-20 A=hilb(4);Y=expm(A)计算结果为：Y=3.2506 1.2068 0.8355 0.6417 1.2068 1.7403 0.5417 0.4288 0.8355 0.5417 1.4100 0.3318 0.6417 0.4288 0.3318 1.2729 函数 log 功能 自然对数，即以 e 为底数的对数。格式 Y=log(X)%对参量 X 中的每一个元素计算自然对数。其中 X 中的元素可以是复数与负数，但由此可能得到意想不到的结果。若 z=x+i*y，则 log 对复数的计算如下：log(z)=log(abs(z)+i*atan2(y,x)例 2-21 下面的语句可以得到无理数的近似值：Pi=abs(log(-1)计算结果为：Pi=3.1416 函数 log10 功能 常用对数，即以 10 为底数的对数。格式 Y=log10(X)%计算 X 中的每一个元素的常用对数，若 X 中出现复数，则可能得到意想不到的结果。例 2-22 MATLAB6.0 数学手册 72L1=log10(realmax)%由此可得特殊变量 realmax 的近似值 L2=log10(eps)%由此可得特殊变量 eps 的近似值 M=magic(4);L3=log10(M)计算结果为：L1=308.2547 L2=-15.6536 L3=1.2041 0.3010 0.4771 1.1139 0.6990 1.0414 1.0000 0.9031 0.9542 0.8451 0.7782 1.0792 0.6021 1.1461 1.1761 0 函数 sort 功能 把输入参量中的元素按从小到大的方向重新排列 格式 B=sort(A)%沿着输入参量 A 的不同维的方向、从小到大重新排列 A 中的元素。A 可以是字符串的、实数的、复数的单元数组。对于 A 中完全相同的元素，则按它们在 A 中的先后位置排列在一块；若 A 为复数的，则按元素幅值的从小到大排列，若有幅值相同的复数元素，则再按它们在区间-,的幅角从小到大排列；若 A 中有元素为NaN，则将它们排到最后。若 A 为向量，则返回从小到大的向量，若 A 为二维矩阵，则按列的方向进行排列；若 A 为多维数组，sort(A)把沿着第一非单元集的元素象向量一样进行处理。B=sort(A,dim)%沿着矩阵 A（向量的、矩阵的或多维的）中指定维数dim 方向重新排列 A 中的元素。B,INDEX=sort(A,)%输出参量 B 的结果如同上面的情形，输出 INDEX 是一等于 size(A)的数组，它的每一列是与 A 中列向量的元素相对应的置换向量。若 A 中有重复出现的相同的值，则返回保存原来相对位置的索引。例 2-23 A=-1.9,-0.2,3.1415926,5.6,7.0,2.4+3.6i;B1,INDEX=sort(A)M=magic(4);B2=sort(M)计算结果为：B1=Columns 1 through 4 -0.2000 -1.9000 3.1416 2.4000+3.6000i Columns 5 through 6 5.6000 7.0000 INDEX=2 1 3 6 4 5 B2=4 2 3 1 5 7 6 8 9 11 10 12 第 2 章 数值计算与数据分析 73 16 14 15 13 函数 abs 功能 数值的绝对值与复数的幅值 格式 Y=abs(X)%返回参量 X 的每一个分量的绝对值；若 X 为复数的，则返回每一分量的幅值：abs(X)=sqrt(real(X).2+imag(X).2)。例 2-24 A=-1.9,-0.2,3.1415926,5.6,7.0,2.4+3.6i;Y=abs(A)计算结果为：Y=1.9000 0.2000 3.1416 5.6000 7.0000 4.3267 函数 conj 功能 复数的共轭值 格式 ZC=conj(Z)%返回参量 Z 的每一个分量的共轭复数：conj(Z)=real(Z)-i*imag(Z)函数 imag 功能 复数的虚数部分 格式 Y=imag(Z)%返回输入参量 Z 的每一个分量的虚数部分。例 2-25 imag(2+3i)计算结果为：ans=3 函数 real 功能 复数的实数部分。格式 Y=real(Z)%返回输入参量 Z 的每一个分量的实数部分。例 2-26 real(2+3i)计算结果为：ans=2 函数 angle 功能 复数的相角 格式 P=angle(Z)%返回输入参量 Z 的每一复数元素的、单位为弧度的相角，其值在区间-，上。说明 angle(z)=imag(log(z)=atan2(imag(z),real(z)例 2-27 Z=1-i,2+i,3-i,4+i;1+2i,2-2i,3+2i,4-2i;1-3i,2+3i,3-3i,4+3i;1+4i,2-4i,3+4i,4-4i;P=angle(Z)计算结果为：P=MATLAB6.0 数学手册 74 -0.7854 0.4636 -0.3218 0.2450 1.1071 -0.7854 0.5880 -0.4636 -1.2490 0.9828 -0.7854 0.6435 1.3258 -1.1071 0.9273 -0.7854 函数 complex 功能 用实数与虚数部分创建复数 格式 c=complex(a,b)%用两个实数 a，b 创建复数 c=a+bi。输出参量 c 与 a、b 同型（同为向量、矩阵、或多维阵列）。该命令比下列形式的复数输入更有用：a+i*b 或 a+j*b 因为 i 和 j 可能被用做其他的变量(不等于 sqrt(-1)，或者 a 和 b 不是双精度的。c=complex(a)%输入参量 a 作为输出复数 c 的实部，其虚部为 0：c=a+0*i。例 2-28 a=uint8(1;2;3;4);b=uint8(4;3;2;1);c=complex(a,b)计算结果为：c=1.0000+4.0000i 2.0000+3.0000i 3.0000+2.0000i 4.0000+1.0000i 函数 mod 功能 模数（带符号的除法余数）用法 M=mod(X,Y)%输入参量 X、Y 应为整数，此时返回余数 X-Y.*floor(X./Y)，若Y0，或者是X。若运算数x与y有相同的符号，则mod(X,Y)等于 rem(X,Y)。总之，对于整数 x,y，有：mod(-x,y)=rem(-x,y)+y。若输入为实数或复数，由于浮点数在计算机上的不精确表示，该操作将导致不可预测的结果。例 2-29 M1=mod(13,5)M2=mod(1:5,3)M3=mod(magic(3),3)计算结果为：M1=3 M2=1 2 0 1 2 M3=2 1 0 0 2 1 1 0 2 函数 nchoosek 功能 二项式系数或所有的组合数。该命令只有对 nC=nchoosek(2:2:10,4)计算结果为：C=2 4 6 8 2 4 6 10 2 4 8 10 2 6 8 10 4 6 8 10 函数 rand 功能 生成元素均匀分布于(0,1)上的数值与阵列 用法 Y=rand(n)%返回 n*n 阶的方阵 Y，其元素均匀分布于区间(0,1)。若 n 不是一标量，在显示一出错信息。Y=rand(m,n)、Y=rand(m n)%返回阶数为 m*n 的，元素均匀分布于区间(0,1)上矩阵 Y。Y=rand(m,n,p,)、Y=rand(m n p)%生成阶数 m*n*p*的，元素服从均匀分布的多维随机阵列 Y。Y=rand(size(A)%生成一与阵列 A 同型的随机均匀阵列 Y rand%该命令在每次单独使用时，都返回一随机数（服从均匀分布）。s=rand(state)%返回一有 35 元素的列向量 s，其中包含均匀分布生成器的当前状态。该改变生成器的当前的状态，见表 2-1。表 2-1 命 令 含 义 Rand(state,s)设置状态为 s Rand(state,0)设置生成器为初始状态 Rand(state,k)设置生成器第 k 个状态(k 为整数)Rand(state,sum(100*clock)设置生成器在每次使用时的状态都不同（因为 clock 每次都不同）例：R1=rand(4,5)a=10;b=50;R2=a+(b-a)*rand(5)%生成元素均匀分布于(10,50)上的矩阵 计算结果可能为：R1=0.6655 0.0563 0.2656 0.5371 0.6797 0.3278 0.4402 0.9293 0.5457 0.6129 0.6325 0.4412 0.9343 0.9394 0.3940 0.5395 0.6501 0.5648 0.7084 0.2206 R2=33.6835 19.8216 36.9436 49.6289 46.4679 18.5164 34.2597 15.3663 31.0549 49.0377 19.0026 37.1006 33.6046 39.5361 13.9336 12.4641 12.9804 35.5420 23.2916 46.8304 28.5238 48.7418 49.0843 13.0512 10.9265 函数 randn MATLAB6.0 数学手册 76功能 生成元素服从正态分布（N(0,1)）的数值与阵列 格式 Y=randn(n)%返回 n*n 阶的方阵 Y，其元素服从正态分布 N(0,1)。若 n 不是一标量，则显示一出错信息。Y=randn(m,n)、Y=randn(m n)%返回阶数为 m*n 的，元素均匀分布于区间(0,1)上矩阵 Y。Y=randn(m,n,p,)、Y=randn(m n p)%生成阶数 m*n*p*的，元素服从正态分布的多维随机阵列 Y。Y=randn(size(A)%生成一与阵列 A 同型的随机正态阵列 Y randn%该命令在每次单独使用时，都返回一随机数（服从正态分布）。s=randn(state)%返回一有 2 元素的向量 s，其中包含正态分布生成器的当前状态。该改变生成器的当前状态，见表 2-2。表 2-2 命 令 含 义 randn(state,s)设置状态为 s randn(state,0)设置生成器为初始状态 rand(state,k)设置生成器第 k 个状态(k 为整数)rand(state,sum(100*clock)设置生成器在每次使用时的状态都不同（因为 clock 每次都不同）例：R1=rand(4,5)R2=0.6+sqrt(0.1)*randn(5)计算结果可能为：R1=0.2778 0.2681 0.5552 0.5167 0.8821 0.2745 0.3710 0.1916 0.3385 0.5823 0.9124 0.5129 0.4164 0.2993 0.0550 0.4125 0.2697 0.1508 0.9370 0.5878 R2=0.4632 0.9766 0.5410 0.6360 0.6931 0.0733 0.9760 0.8295 0.9373 0.1775 0.6396 0.5881 0.4140 0.6187 0.8259 0.6910 0.7035 1.2904 0.5698 1.1134 0.2375 0.6552 0.5569 0.3368 0.3812 2.2 插值、拟合与查表插值、拟合与查表 插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量 x 与因变量 y 的函数 y=f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。如何根据观测点的值，构造一个比较简单的函数 y=(x)，使函数在观测点的值等于已知的数值或导数值。用简单函数 y=(x)在点 x 处的值来估计未知函数 y=f(x)在 x 点的值。寻找这样的函数(x)，办法是很多的。(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式；(x)可以是任意光滑（任意阶导数连续）的函数或是分段函数。函数类的不同，自然地有不同的逼近效果。在许多应用中，通常要用一个解析函数（一、二元函数）来描述观测数据。第 2 章 数值计算与数据分析 77根据测量数据的类型：1测量值是准确的，没有误差。2测量值与真实值有误差。这时对应地有两种处理观测数据方法：1插值或曲线拟合。2回归分析（假定数据测量是精确时，一般用插值法，否则用曲线拟合）。MATLAB 中提供了众多的数据处理命令。有插值命令，有拟合命令，有查表命令。2.2.1 插值命令插值命令 命令 1 interp1 功能 一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数 f(x)由所给数据决定。各个参量之间的关系示意图为图 2-14。f(x)x：原始数据点Y：原始数据点xi：插值点Yi：插值点 图 2-14 数据点与插值点关系示意图 格式 yi=interp1(x,Y,xi)%返回插值向量 yi，每一元素对应于参量 xi，同时由向量 x 与 Y 的内插值决定。参量 x 指定数据 Y 的点。若 Y 为一矩阵，则按 Y 的每列计算。yi 是阶数为length(xi)*size(Y,2)的输出矩阵。yi=interp1(Y,xi)%假定 x=1:N，其中 N 为向量 Y 的长度，或者为矩阵Y 的行数。yi=interp1(x,Y,xi,method)%用指定的算法计算插值：nearest：最近邻点插值，直接完成计算；linear：线性插值（缺省方式），直接完成计算；spline：三次样条函数插值。对于该方法，命令 interp1 调用函数 spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函数。命令 spline 用它们执行三次样条函数插值；pchip：分段三次 Hermite 插值。对于该方法，命令 interp1 调用函数 pchip，用于对向量 x 与 y 执行分段三次内插值。该方法保留单调性与数据的外形；cubic：与pchip操作相同；v5cubic：在 MATLAB 5.0 中的三次插值。对于超出 x 范围的 xi 的分量，使用方法nearest、linear、v5cubic的插值算法，相应地将返回 NaN。对其他的方法，interp1 将对超出的分量执行外插值算法。MATLAB6.0 数学手册 78yi=interp1(x,Y,xi,method,extrap)%对于超出 x 范围的 xi 中的分量将执行特殊的外插值法 extrap。yi=interp1(x,Y,xi,method,extrapval)%确定超出 x 范围的 xi 中的分量的外插值extrapval，其值通常取 NaN 或 0。例 2-31 x=0:10;y=x.*sin(x);xx=0:.25:10;yy=interp1(x,y,xx);plot(x,y,kd,xx,yy)插值图形为图 2-15。例 2-32 year=1900:10:2010;product=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 256.344 267.893;p1995=interp1(year,product,1995)x=1900:1:2010;y=interp1(year,product,x,pchip);plot(year,product,o,x,y)插值结果为：p1995=252.9885 插值图形为图 2-16。图 2-15 一元函数插值图形 图 2-16 离散数据的一维插值图 命令 2 interp2 功能 二维数据内插值（表格查找）格式 ZI=interp2(X,Y,Z,XI,YI)%返回矩阵 ZI，其元素包含对应于参量 XI 与 YI（可以是向量、或同型矩阵）的元素，即 Zi(i,j)Xi(i,j),yi(i,j)。用户可以输入行向量和列向量 Xi 与Yi，此时，输出向量 Zi 与矩阵 meshgrid(xi,yi)是同型的。同时取决于由输入矩阵 X、Y 与 Z 确定的二维函数 Z=f(X,Y)。参量 X 与 Y 必须是单调的，且相同的划分格式，就像由命令 meshgrid 生成的一样。若 Xi与Yi中有在X与Y范围之外的点，则相应地返回nan（Not a Number）。ZI=interp2(Z,XI,YI)%缺省地，X=1:n、Y=1:m，其中m,n=size(Z)。再按第 2 章 数值计算与数据分析 79第一种情形进行计算。ZI=interp2(Z,n)%作 n 次递归计算，在 Z 的每两个元素之间插入它们的二维插值，这样，Z 的阶数将不断增加。interp2(Z)等价于 interp2(z,1)。ZI=interp2(X,Y,Z,XI,YI,method)%用指定的算法 method 计算二维插值：linear：双线性插值算法（缺省算法）；nearest：最临近插值；spline：三次样条插值；cubic：双三次插值。例 2-33：X,Y=meshgrid(-3:.25:3);Z=peaks(X,Y);XI,YI=meshgrid(-3:.125:3);ZZ=interp2(X,Y,Z,XI,YI);surfl(X,Y,Z);hold on;surfl(XI,YI,ZZ+15)axis(-3 3-3 3-5 20);shading flat hold off 插值图形为图 2-17。例 2-34 years=1950:10:1990;service=10:10:30;wage=150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243;w=interp2(service,years,wage,15,1975)插值结果为：w=190.6288 命令 3 interp3 功能 三维数据插值（查表）格式 VI=interp3(X,Y,Z,V,XI,YI,ZI)%找出由参量X,Y,Z决定的三元函数V=V(X,Y,Z)在点（XI,YI,ZI）的值。参量 XI,YI,ZI 是同型阵列或向量。若向量参量 XI,YI,ZI 是不同长度，不同方向（行或列）的向量，这时输出参量 VI 与 Y1,Y2,Y3 为同型矩阵。其中 Y1,Y2,Y3 为用命令meshgrid(XI,YI,ZI)生成的同型阵列。若插值点(XI,YI,ZI)中有位于点(X,Y,Z)之外的点，则相应地返回特殊变量值 NaN。VI=interp3(V,XI,YI,ZI)%缺省地，X=1:N，Y=1:M，Z=1:P，其中，M,N,P=size(V)，再按上面的情形计算。VI=interp3(V,n)%作 n 次递归计算，在 V 的每两个元素之间插入它们的三维插值。这样，V 的阶数将不断增加。interp3(V)等价于 interp3(V,1)。图 2-17 二维插值图 MATLAB6.0 数学手册 80VI=interp3(,method)%用指定的算法 method 作插值计算：linear：线性插值（缺省算法）；cubic：三次插值；spline：三次样条插值；nearest：最邻近插值。说明 在所有的算法中，都要求 X,Y,Z 是单调且有相同的格点形式。当 X,Y,Z 是等距且单调时，用算法*linear，*cubic，*nearest，可得到快速插值。例 2-35 x,y,z,v=flow(20);xx,yy,zz=meshgrid(.1:.25:10,-3:.25:3,-3:.25:3);vv=interp3(x,y,z,v,xx,yy,zz);slice(xx,yy,zz,vv,6 9.5,1 2,-2.2);shading interp;colormap cool 插值图形为图 2-18。图 2-18 三维插值图 命令 4 interpft 功能 用快速 Fourier 算法作一维插值 格式 y=interpft(x,n)%返回包含周期函数 x 在重采样的 n 个等距的点的插值 y。若length(x)=m，且 x 有采样间隔 dx，则新的 y 的采样间隔dy=dx*m/n。注意的是必须 nm。若 x 为一矩阵，则按 x 的列进行计算。返回的矩阵 y 有与 x 相同的列数，但有 n 行。y=interpft(x,n,dim)%沿着指定的方向 dim 进行计算 命令 5 griddata 功能 数据格点 格式 ZI=griddata(x,y,z,XI,YI)%用二元函数 z=f(x,y)的曲面拟合有不规则的数据向量 x,y,z。griddata 将返回曲面 z 在点（XI,YI）处的插值。曲面总是经过这些数据点（x,y,z）的。输入参量（XI,YI）通常