1、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),
2、,2e)zsin(izi2e)sn(z例 2-1x = -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 中每一个
3、元素的反双曲正弦函数值 Y第 2 章 数值计算与数据分析63说明 反正弦函数与反双曲正弦函数的定义为: ,)z1iln(zsia2)z1ln(zsiha2例 2-2x = -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)并不是精确的零,而是与浮点精度有关的无穷小
4、量eps,因为 pi 仅仅是精确值 浮点近似的表示值而已。Y = sinh(X) %计算参量 X 的双曲余弦值 Y说明 若 X 为复数 z= x+iy,则函数定义为:cos(x+iy) = cos(x)*cos(y) + i*sin(x)*sin(y),2ezcosizi2ecoshz例 2-3x = -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) %返回参
5、量 X(可以是向量、矩阵)中每一个元素的反余弦函数值 Y。若 X 中有的分量处于-1,1之间,则 Y = acos(X)对应的分量处于0, 之间,若 X 中有分量在区间-1,1 之外,则 Y = acos(X)对应的分量为复数。Y = asinh(X) %返回参量 X 中每一个元素的反双曲余弦函数 Y说明 反余弦函数与反双曲余弦函数定义为: ,)z1iln(izcosa2)1zln(zcosha2例 2-4x = -1:.01:1; plot(x,acos(x)x = -5:.01:5; plot(x,acosh(x)图形结果为图 2-4。图 2-4 反余弦函数与反双曲余弦函数图函数 tan、
6、tanh功能 正切函数与双曲正切函数格式 Y = tan(X) %计算参量 X(可以是向量、矩阵,元素可以是复数)中每一个角度分量的正切值 Y,所有角度分量的单位为弧度。我们要指出的是,tan(pi/2)并不是精确的零,而是与浮点精度有关的无穷小量eps,因为 pi 仅仅是精确值 浮点近似的表示值而已。Y = tanh(X) %返回参量 X 中每一个元素的双曲正切函数值 Y例 2-5x = (-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 正切函数
7、与双曲正切函数图第 2 章 数值计算与数据分析65函数 atan、atanh功能 反正切函数与反双曲正切函数格式 Y = atan(X) %返回参量 X(可以是向量、矩阵)中每一个元素的反正切函数值 Y。若 X 中有的分量为实数,则 Y = atan(X)对应的分量处于-/2,/2之间。Y = atanh(X) %返回参量 X 中每一个元素的反双曲正切函数值 Y。说明 反正切函数与反双曲正切函数定义为: ,ziln2taz1ln2tah例 2-6x = -20:0.01:20; plot(x,atan(x)x = -0.99:0.01:0.99; plot(x,atanh(x)图形结果为图 2
8、-6。图 2-6 反正切函数与反双曲正切函数图函数 cot、coth功能 余切函数与双曲余切函数格式 Y = cot(X) %计算参量 X(可以是向量、矩阵,元素可以是复数)中每一个角度分量的余切值 Y,所有角度分量的单位为弧度。Y = coth(X) %返回参量 X 中每一个元素的双曲余切函数值 Y例 2-7x1 = -pi+0.01:0.01:-0.01; % 去掉奇点 x = 0x2 = 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 数
9、学手册66图 2-7 余切函数与双曲余切函数图函数 acot、acoth功能 反余切函数与反双曲余切函数格式 Y = acot(X) %返回参量 X(可以是向量、矩阵)中每一个元素的反余切函数YY = acoth(X) %返回参量 X 中每一个元素的反双曲余切函数值 Y例 2-8x1 = -2*pi:pi/30:-0.1; x2 = 0.1:pi/30:2*pi; % 去掉奇异点 x = 0plot(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
10、。图 2-8 反余切函数与反双曲余切函数图函数 sec、sech功能 正割函数与双曲正割函数格式 Y = sec(X) %计算参量 X(可以是向量、矩阵,元素可以是复数)中每一个角度分量的正割函数值 Y,所有角度分量的单位为弧度。我们要指出的是,sec(pi/2)并不是无穷大,而是与浮点精度有关的无穷小量 eps 的倒数,因为 pi 仅仅是精确值 浮点近似的表示值而已。Y = sech(X) %返回参量 X 中每一个元素的双曲正割函数值 Y例 2-9x1 = -pi/2+0.01:0.01:pi/2-0.01; % 去掉奇异点 x = pi/2x2 = pi/2+0.01:0.01:(3*pi
11、/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(可以是向量、矩阵)中每一个元素的反正割函数值 YY = asech(X) %返回参量 X 中每一个元素的反双曲正割函数值 Y例 2-10x1 = -5:0.01:-1; x2 = 1:0.01:5; plot(x1,asec(x1),x2,asec(x2)x = 0
12、.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-11x1 = -pi+0.01:0.01:-0.01; x2 = 0.01:0.01:pi-0.01; % 去掉奇异点 x=0plot(x1,csc(x1),x2,csc(x2)plot(x1,csch(x1)
13、,x2,csch(x2)图形结果为图 2-11。MATLAB6.0 数学手册68图 2-11 余割函数与双曲余割函数图函数 acsc、acsch功能 反余割函数与反双曲余割函数。格式 Y = asec(X) %返回参量 X(可以是向量、矩阵)中每一个元素的反余割函数值 YY = asech(X) %返回参量 X 中每一个元素的反双曲余割函数值 Y例 2-12x1 = -10:0.01:-1.01; x2 = 1.01:0.01:10; % 去掉奇异点 x = 1plot(x1,acsc(x1),x2,acsc(x2)x1 = -20:0.01:-1; x2 = 1:0.01:20; plot(
14、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-13z=1+2i;r = abs(z);theta = atan2(imag(z),real(z) z = r *exp(i *theta)feathe
15、r(z);hold on第 2 章 数值计算与数据分析69t=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.1071z =1.0000 + 2.0000i图形结果为图 2-13。图 2-13 四象限的反正切函数图2.1.2 其他常用函数函数 fix功能 朝零方向取整格式 B = fix(A) %对 A 的每一个元素朝零的方向取整数部分,返回与 A 同维的数组。对于复数参量 A,则返回一复数,其分量的实数与虚数部分分别取原复数的、朝零方向的整数部分。例 2-1
16、4A = -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-15A = -1.9, -0.2, 3.1415926
17、, 5.6, 7.0, 2.4+3.6i;Y = round(A)MATLAB6.0 数学手册70计算结果为: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-16A = -1.9, -0.2, 3.1415926, 5.
18、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)返回的结果与
19、 mod(X,Y)相同,不然,若 X 为正数,则 rem(-X,Y) = mod(-X,Y) - Y。该命令返回的结果在区间0 ,sign(X)*abs(Y),若 Y 中有零分量,则相应地返回 NaN。例 2-17X = 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-18A = -1.9, -
20、0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i;B = ceil(A)计算结果为:第 2 章 数值计算与数据分析71B =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-19A = -1
21、.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.2704Columns 5 through 6 1.0966 -0.0099 - 0.0049i函数 expm功能 求矩阵的以 e 为底数的指数函数格式 Y = expm(X) %计算以 e 为底数、x 的每一个元素为指数的指数函数值。若矩阵 x 有小于等于零的特征值,则返回复数的结果。说明 该函数为一内建函数,它有三种计算算法:(1)使用文件 expm1.m 中的用比例法
22、与二次幂算法得到的 Pad 近似值;(2)使用 Taylor 级数近似展开式计算,这种计算在文件 expm2.m 中。但这种一般计算方法是不可取的,通常计算是缓慢且不精确的;(3)在文件 expm3.m 中,先是将矩阵对角线化,再把函数计算出相应的的特征向量,最后转换过来。但当输入的矩阵没有与矩阵阶数相同的特征向量个数时,就会出现错误。例 2-20A=hilb(4);Y = expm(A)计算结果为:Y =3.2506 1.2068 0.8355 0.64171.2068 1.7403 0.5417 0.42880.8355 0.5417 1.4100 0.33180.6417 0.4288
23、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 中的每
24、一个元素的常用对数,若 X 中出现复数,则可能得到意想不到的结果。例 2-22L1 = log10(realmax) % 由此可得特殊变量 realmax 的近似值L2 = log10(eps) % 由此可得特殊变量 eps 的近似值M = magic(4);L3 = log10(M)计算结果为:L1 =308.2547L2 =-15.6536L3 =1.2041 0.3010 0.4771 1.11390.6990 1.0414 1.0000 0.90310.9542 0.8451 0.7782 1.07920.6021 1.1461 1.1761 0函数 sort功能 把输入参量中的元素按
25、从小到大的方向重新排列格式 B = sort(A) %沿着输入参量 A 的不同维的方向、从小到大重新排列 A 中的元素。A 可以是字符串的、实数的、复数的单元数组。对于 A 中完全相同的元素,则按它们在 A 中的先后位置排列在一块;若 A 为复数的,则按元素幅值的从小到大排列,若有幅值相同的复数元素,则再按它们在区间-, 的幅角从小到大排列;若 A 中有元素为 NaN,则将它们排到最后。若 A 为向量,则返回从小到大的向量,若 A 为二维矩阵,则按列的方向进行排列;若 A 为多维数组,sort(A)把沿着第一非单元集的元素象向量一样进行处理。B = sort(A,dim) %沿着矩阵 A(向量
26、的、矩阵的或多维的)中指定维数dim 方向重新排列 A 中的元素。B,INDEX = sort(A,) %输出参量 B 的结果如同上面的情形,输出 INDEX是一等于 size(A)的数组,它的每一列是与 A 中列向量的元素相对应的置换向量。若 A 中有重复出现的相同的值,则返回保存原来相对位置的索引。例 2-23A = -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.200
27、0 -1.9000 3.1416 2.4000 + 3.6000iColumns 5 through 6 5.6000 7.0000INDEX =2 1 3 6 4 5B2 =4 2 3 15 7 6 89 11 10 1216 14 15 13函数 abs功能 数值的绝对值与复数的幅值格式 Y = abs(X) %返回参量 X 的每一个分量的绝对值;若 X 为复数的,则返回每一分量的幅值:abs(X) = sqrt(real(X).2+imag(X).2) 。例 2-24A = -1.9, -0.2, 3.1415926, 5.6, 7.0, 2.4+3.6i;Y = abs(A)计算结果为
28、: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-25imag(2+3i)计算结果为:ans = 3函数 real功能 复数的实数部分。格式 Y = real(Z) %返回输入参量 Z 的每一个分量的实数部分。例 2-26real(2+3i)计算结果为:ans =2函数 angl
29、eMATLAB6.0 数学手册74功能 复数的相角格式 P = angle(Z) %返回输入参量 Z 的每一复数元素的、单位为弧度的相角,其值在区间-,上。说明 angle(z) = imag (log(z) = atan2 (imag(z),real(z)例 2-27Z =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.24501.1071 -0.7854 0.5880 -0.4636-1.2
30、490 0.9828 -0.7854 0.64351.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-28a
31、 = uint8(1;2;3;4);b = uint8(4;3;2;1);c = complex(a,b)计算结果为:c =1.0000 + 4.0000i2.0000 + 3.0000i3.0000 + 2.0000i4.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。若输入为实数或复数
32、,由于浮点数在计算机上的不精确表示,该操作将导致不可预测的结果。例 2-29M1 = mod(13,5) M2 = mod(1:5,3) M3 = mod(magic(3),3)第 2 章 数值计算与数据分析75计算结果为:M1 =3M2 =1 2 0 1 2M3 =2 1 00 2 11 0 2函数 nchoosek功能 二项式系数或所有的组合数。该命令只有对 nC = nchoosek(2:2:10,4)计算结果为:C =2 4 6 82 4 6 102 4 8 102 6 8 104 6 8 10函数 rand功能 生成元素均匀分布于(0,1)上的数值与阵列用法 Y = rand(n)
33、%返回 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 同型的随机均匀阵列 Yrand %该命令在每次单独使用时,都返回一随机数(服从均匀分布)。s = rand(state) %返回一有 35 元素的列向量 s,其中包含均匀分布生成器的当前状态
34、。该改变生成器的当前的状态,见表 2-1。表 2-1命 令 含 义Rand(state,s) 设置状态为 sRand(state,0) 设置生成器为初始状态Rand(state,k) 设置生成器第 k 个状态(k 为整数)Rand(state,sum(100*clock) 设置生成器在每次使用时的状态都不同(因为 clock 每次都不同)MATLAB6.0 数学手册76例:R1 = rand(4,5)a = 10; b = 50;R2 = a + (b-a) * rand(5) % 生成元素均匀分布于 (10,50)上的矩阵计算结果可能为:R1 =0.6655 0.0563 0.2656 0.
35、5371 0.67970.3278 0.4402 0.9293 0.5457 0.61290.6325 0.4412 0.9343 0.9394 0.39400.5395 0.6501 0.5648 0.7084 0.2206R2 =33.6835 19.8216 36.9436 49.6289 46.467918.5164 34.2597 15.3663 31.0549 49.037719.0026 37.1006 33.6046 39.5361 13.933612.4641 12.9804 35.5420 23.2916 46.830428.5238 48.7418 49.0843 13.
36、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 同型的随机正态阵列 Yrandn %该命令在每次单独使用
37、时,都返回一随机数(服从正态分布) 。s = randn(state) %返回一有 2 元素的向量 s,其中包含正态分布生成器的当前状态。该改变生成器的当前状态,见表 2-2。表 2-2命 令 含 义randn(state,s) 设置状态为 srandn(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.
38、2778 0.2681 0.5552 0.5167 0.88210.2745 0.3710 0.1916 0.3385 0.58230.9124 0.5129 0.4164 0.2993 0.05500.4125 0.2697 0.1508 0.9370 0.5878R2 =0.4632 0.9766 0.5410 0.6360 0.69310.0733 0.9760 0.8295 0.9373 0.17750.6396 0.5881 0.4140 0.6187 0.82590.6910 0.7035 1.2904 0.5698 1.11340.2375 0.6552 0.5569 0.336
39、8 0.3812第 2 章 数值计算与数据分析772.2 插值、拟合与查表插值法是实用的数值方法,是函数逼近的重要方法。在生产和科学实验中,自变量 x与因变量 y 的函数 y = f(x)的关系式有时不能直接写出表达式,而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时,需要估计函数值在该点的值。如何根据观测点的值,构造一个比较简单的函数 y=(x),使函数在观测点的值等于已知的数值或导数值。用简单函数 y=(x)在点 x 处的值来估计未知函数 y=f(x)在 x 点的值。寻找这样的函数 (x),办法是很多的。(x)可以是一个代数多项式,或是三角多项式,也可以是有理分式;
40、(x)可以是任意光滑(任意阶导数连续)的函数或是分段函数。函数类的不同,自然地有不同的逼近效果。在许多应用中,通常要用一个解析函数(一、二元函数)来描述观测数据。根据测量数据的类型:1测量值是准确的,没有误差。2测量值与真实值有误差。这时对应地有两种处理观测数据方法:1插值或曲线拟合。2回归分析(假定数据测量是精确时,一般用插值法,否则用曲线拟合) 。MATLAB 中提供了众多的数据处理命令。有插值命令,有拟合命令,有查表命令。2.2.1 插值命令命令 1 interp1功能 一维数据插值(表格查找) 。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数 f(x)由所
41、给数据决定。各个参量之间的关系示意图为图 2-14。 f(x) x: 原 始 数 据 点 Y: 原 始 数 据 点 i: 插 值 点 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 的行数。MATLAB6.0 数学
42、手册78yi = interp1(x,Y,xi,method) %用指定的算法计算插值:nearest:最近邻点插值,直接完成计算;linear:线性插值(缺省方式) ,直接完成计算;spline:三次样条函数插值。对于该方法,命令 interp1 调用函数spline、ppval、mkpp、umkpp 。这些命令生成一系列用于分段多项式操作的函数。命令 spline 用它们执行三次样条函数插值;pchip:分段三次 Hermite 插值。对于该方法,命令 interp1 调用函数 pchip,用于对向量 x 与 y 执行分段三次内插值。该方法保留单调性与数据的外形;cubic:与 pchip
43、操作相同;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-31x = 0:10;
44、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);p
45、lot(year,product,o,x,y)插值结果为:p1995 =252.9885插值图形为图 2-16。第 2 章 数值计算与数据分析79图 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、
46、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
47、,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 flathold off插值图形为图 2-17。例 2-34 图 2-17 二维插值图MATLAB6.0 数学手册80years = 1950:10:1990;service = 10:10:30;wage = 150.697 199.592 187.625179.323 195.072 250.287203.212 179.092 322.767226.505 153.706 426.730249.633 120.281 598.243;w = interp2(service,years,wage,15,1975)插值结果为:w =190.6288命令 3 interp3功能 三维数据插值(查表)格式 VI = interp3(X,