MATLAB数学手册教程_第2章_数值计算与数据分析.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=第 2 章 数值计算与数据分析 63asin(X)对应的分量为复数。Y=asinh(X)%返回参量 X 中每一个元素的反双曲正弦函数值 Y 说明 反正弦函数与反双曲正弦函数的定义为:)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 余弦函数与双曲余弦函数图 MATLAB6.0 数学手册 64函数 acos、acosh 功能 反余弦函数与反双曲余弦函数 格式 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 章 数值计算与数据分析 65 图 2-5 正切函数与双曲正切函数图 函数 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。MATLAB6.0 数学手册 66 图 2-7 余切函数与双曲余切函数图 函数 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 章 数值计算与数据分析 67 图 2-9 正割函数与双曲正割函数图 函数 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。MATLAB6.0 数学手册 68 图 2-11 余割函数与双曲余割函数图 函数 acsc、acsch 功能 反余割函数与反双曲余割函数。格式 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 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;MATLAB6.0 数学手册 70Y=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,则返回一复数,其分量的实数与虚数部分分别取原复数的、朝负无穷大方向的整数部分。例 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)第 2 章 数值计算与数据分析 71计算结果为: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-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)MATLAB6.0 数学手册 72例 2-21 下面的语句可以得到无理数的近似值:Pi=abs(log(-1)计算结果为:Pi=3.1416 函数 log10 功能 常用对数,即以 10 为底数的对数。格式 Y=log10(X)%计算 X 中的每一个元素的常用对数,若 X 中出现复数,则可能得到意想不到的结果。例 2-22 L1=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);第 2 章 数值计算与数据分析 73B2=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 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 MATLAB6.0 数学手册 74功能 复数的相角 格式 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=-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)计算结果为:第 2 章 数值计算与数据分析 75M1=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;MATLAB6.0 数学手册 76R2=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 功能 生成元素服从正态分布(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 插值、拟合与查表插值、拟合与查表 第 2 章 数值计算与数据分析 77插值法是实用的数值方法,是函数逼近的重要方法。在生产和科学实验中,自变量 x 与因变量 y 的函数 y=f(x)的关系式有时不能直接写出表达式,而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时,需要估计函数值在该点的值。如何根据观测点的值,构造一个比较简单的函数 y=(x),使函数在观测点的值等于已知的数值或导数值。用简单函数 y=(x)在点 x 处的值来估计未知函数 y=f(x)在 x 点的值。寻找这样的函数(x),办法是很多的。(x)可以是一个代数多项式,或是三角多项式,也可以是有理分式;(x)可以是任意光滑(任意阶导数连续)的函数或是分段函数。函数类的不同,自然地有不同的逼近效果。在许多应用中,通常要用一个解析函数(一、二元函数)来描述观测数据。根据测量数据的类型: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、MATLAB6.0 数学手册 78mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函数。命令 spline 用它们执行三次样条函数插值;pchip:分段三次 Hermite 插值。对于该方法,命令 interp1 调用函数 pchip,用于对向量 x 与 y 执行分段三次内插值。该方法保留单调性与数据的外形;cubic:与pchip操作相同;v5cubic:在 MATLAB 5.0 中的三次插值。对于超出 x 范围的 xi 的分量,使用方法nearest、linear、v5cubic的插值算法,相应地将返回 NaN。对其他的方法,interp1 将对超出的分量执行外插值算法。yi=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(可第 2 章 数值计算与数据分析 79以是向量、或同型矩阵)的元素,即 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)。再按第一种情形进行计算。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 是同型阵列或向量。若向量图 2-17 二维插值图 MATLAB6.0 数学手册 80参量 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)。VI=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 第 2 章 数值计算与数据分析 81功能 数据格点 格式 ZI=griddata(x,y,z,XI,YI)%用二元函数 z=f(x,y)的曲面拟合有不规则的数据向量 x,y,z。griddata 将返回曲面 z 在点(XI,YI)处的插值。曲面总是经过这些数据点(x,y,z)的。输入参量(XI,YI)通常是规则的格点(像用命令 meshgrid 生成的一样)。XI 可以是一行向量,这时 XI 指定一有常数列向量的矩阵。类似地,YI 可以是一列向量,它指定一有常数行向量的矩阵。XI