matlab第2 数值数组 (2).ppt

第二讲 数值数组及其运算 * 1 一、一维数组的创建和寻访 1、一维数组的创建 （1）逐个元素输入法 例：x=2 pi/2 sqrt(3) 3+5i x= 2.0000 1.5708 1.7321 3.0000+5.0000i （2）冒号生成法 该法通过“步长”设定，生成一维“行”数组的方法。 例：x=a:inc:b %inc是采样点的步长 （3）定数线性采样法 该法在设定的“总点数”下，均匀采样生成一维“行”数组。 例：x=linspace(a,b,n) %a,b分别是生成数组的第一和最 后一个元素，n是采样总点数。 * 2 （4）定数对数采样法 该法在设定的“总点数”下，经“常用对数”采样生成一维 “行”数组 例：x=logspace(a,b,n) %a,b分别代表生成数组的第一 和最后元素分别为：10a，10b，n是采样点数。 2、一维数组的子数组寻访和赋值 【例】子数组的寻访。 rand(state,0) %把均匀分布伪随机发生器置为0 x=rand(1,5) %产生（15）的均匀分布随机数组 x = 0.9501 0.2311 0.6068 0.4860 0.8913 x(3)%寻访数组x的第三个 ans = 0.6068 * 3 x(1 2 5) %寻访数组的第一、二、五元素组成的子数组 ans = 0.9501 0.2311 0.8913 x(1:3) %寻访数组前三个元素组成的子数组 ans = 0.9501 0.2311 0.6068 x(3:end) %寻访除三个元素外的全部其它元素。end是最 后一个元素的下标 ans = 0.6068 0.4860 0.8913 x(3:-1:1) %由前3个元素倒排构成的子数组 ans = 0.6068 0.2311 0.9501 x(find(x0.5) %由大于0.5 的元素构成的子数组 ans = 0.9501 0.6068 0.8913 * 4 x(1 2 3 4 4 3 2 1) %对元素可以重复寻访，使所得数组长 度允许大于原数组 ans = Columns 1 through 7 0.9501 0.2311 0.6068 0.4860 0.4860 0.6068 0.2311 Column 8 0.9501 二、二维数组的创建 1、直接输入法 【例】在MATLAB环境下，用下面三条指令创建二维数组C。 a=2.7358; b=33/79; C=1,2*a+i*b,b*sqrt(a);sin(pi/4),a+5*b,3.5+i C = 1.0000 5.4716 + 0.4177i 0.6909 0.7071 4.8244 3.5000 + 1.0000i * 5 注意事项： 整个数组必须以方括号为其首尾 数组的行与行之间必须用分号或回车隔离； 数组元素必须由逗号或空格分隔。 2、利用M文件创建和保存数组 （1）打开文件编辑调试器，并在空白填写框中输入以下内容 % MyMatrix.mCreation and preservation of matrix AM AM=101,102,103,104,105,106,107,108,109;. 201,202,203,204,205,206,207,208,209;. 301,302,303,304,305,306,307,308,309; （2）保存此文件，并文件起名为MyMatrix.m （3）只要在指令窗中运行此文件，数组AM就会自动生成于 Matlab内存中. * 6 3利用冒号表达式建立一个向量 冒号表达式可以产生一个行向量，一般格式是： e1:e2:e3 其中e1为初始值，e2为步长，e3为终止值。 4建立大矩阵 大矩阵可由方括号中的小矩阵或向量建立起来。例如： A=1,2,3;4,5,6;7,8,9; C=A,eye(size(A);ones(size(A),A C = 1 2 3 1 0 0 4 5 6 0 1 0 7 8 9 0 0 1 1 1 1 1 2 3 1 1 1 4 5 6 1 1 1 7 8 9 * 7 三、二维数组元素的标识 1、“全下标”标识 全下标标识由两个下标组成：行下标、列下标 2、“单下标”标识 单下标标识就是由一个下标来指明元素在数组中的位置 3、“逻辑1”标识 寻找数组中所有大于某值的元素的问题。 【例】找出数组 中所有绝对值大于3的元素。 A=zeros(2,5);%预生成一个（25）全零数组 A(:)=-4:5%运用“全元素”法向A赋值 L=abs(A)3%产生与A同维的“01”逻辑值数组 islogical(L)%判断L是否逻辑值数组。输出若为1，则 是 X=A(L)%把L中逻辑值1对应的A元素取出 * 8 A = -4 -2 0 2 4 -3 -1 1 3 5 L = 1 0 0 0 1 0 0 0 0 1 ans = 1 X = -4 4 5 * 9 四、二维数组的子数组寻访和赋值 子数组的寻访和赋值使用说明 A(r,c)它由A的“r指定行”和“c指定列”上的元素组成 A(r,:)它由A的“r指定行”和“全部列”上的元素组成 A(:,c)它由A的“全部行”和“c指定列”上的元素组成 A(:)“单下标全元素”寻访，它由A的各列按自左至右的次序，首 尾相接而生成“一维长列”数组 A(s)“单下标”寻访，生成“s指定”一维数组，s若是“行数组” （或列数组），则A(s)是长度相同的“行数组”（或列数组 A(L)“逻辑1”寻访，生成一维列数组，由与A同样大小的”逻辑数 组L中的“1”元素选出A的对应元素，按单下标次序排成长列 A(r,c)=Sa以“双下标”方式，对子数组A(r,c)进行赋值，Sa的“行宽，列 长”必须与A(r,c)的“行宽、列长”相同 A(:)=D(:)全元素赋值方式，结果保持A的行宽、列长不变，条件是A、D 两个数组的总元素数相等，但行宽、列长不一定相同 A(s)=Sa按单下标方式对A的部分元素重新赋值，结果保持A的行宽、列 长不变，s单下标数组的长度必须与“一维数组”Sa的长度相等 。 * 10 1矩阵元素的提取 通过下标引用矩阵的元素，例如，A(3,2)=8 利用矩阵的序号来引用矩阵元素。矩阵元素的序号就是相 应元素在内存中的排列顺序。 矩阵元素按列存储，先第一列，再第二列，依次类推 。 【例】 A=1,2,3;4,5,6; A(3) ans = 2 序号(Index)与下标(Subscript )是一一对应的。 * 11 2矩阵的拆分 (1)利用冒号表达式获得子矩阵 A(:,j)表示取A矩阵的第j列全部元素； A(i,:)表示A矩阵第i行的全部元素； A(i,j)表示取A矩阵第i行、第j列的元素。 (2) A(i:i+m,:)表示取A矩阵第ii+m行的全部元素； A(:,k:k+m)表示取A矩阵第kk+m列的全部元素 ， A(i:i+m,k:k+m)表示取A矩阵第ii+m行内 ，并在第kk+m列中的所有元素。 (3) 利用一般向量和end运算符来表示矩阵下标，从 而获得子矩阵。end表示某一维的末尾元素下标。 * 12 C(1,4,3:end) ans = 3 1 0 0 1 1 2 3 (4) 利用空矩阵删除矩阵的元素 在MATLAB中，定义 为空矩阵。给变量X赋 空矩阵的语句为X= 。注意，X= 与clear X不 同，clear是将X从工作空间中删除，而空矩阵则 存在于工作空间中，只是维数为0。 C = 1 2 3 1 0 0 4 5 6 0 1 0 7 8 9 0 0 1 1 1 1 1 2 3 1 1 1 4 5 6 1 1 1 7 8 9 * 13 【例】不同赋值方式示例。 A=zeros(2,4) %创建(2*4)的全零数组 A = 0 0 0 0 0 0 0 0 A(:)=1:8 %全元素赋值方式 A = 1 3 5 7 2 4 6 8 s=2 3 5 %产生单下标数组行数组 s = 2 3 5 A(s)%由“单下标行数组”寻访产生A元素组成的行数组 ans = 2 3 5 * 14 Sa=10 20 30 %Sa是长度为3的“列数组” Sa = 10 20 30 A(s)=Sa %单下标方式赋值 A = 1 20 30 7 10 4 6 8 A(:,2 3)=ones(2) %双下标赋值方式：把A的第2、3列元素 全赋为1 A = 1 1 1 7 10 1 1 8 * 15 【例】演示pow2的数组运算性质。 A=1:4;5:8 %生成A(2*4)数组 A = 1 2 3 4 5 6 7 8 pow2(A) %计算2的幂也生成数组 ans = 2 4 8 16 32 64 128 256 五、执行数组运算的常用函数 * 16 六、数组运算和矩阵运算 【例 】两种不同转置的比较 clear;A=zeros(2,3); A(:)=1:6; %全元素赋值 A=A*(1+i) %运用标量与数组乘产生复数矩阵 A_A=A. %数组转置，即非共轭转置 A_M=A %矩阵转置，即共轭转置 A = 1.0000 + 1.0000i 3.0000 + 3.0000i 5.0000 + 5.0000i 2.0000 + 2.0000i 4.0000 + 4.0000i 6.0000 + 6.0000i A_A = 1.0000 + 1.0000i 2.0000 + 2.0000i 3.0000 + 3.0000i 4.0000 + 4.0000i 5.0000 + 5.0000i 6.0000 + 6.0000i A_M = 1.0000 - 1.0000i 2.0000 - 2.0000i 3.0000 - 3.0000i 4.0000 - 4.0000i 5.0000 - 5.0000i 6.0000 - 6.0000i * 17 指令含义指令含义 A.非共轭转置，相当于conj(A)A-B对应元素相减 A=s把标量s赋给A的每个元素A.B对应元素相乘 s+B标量s分别与B元素之和A./BA的元素被B的对应元素除 s-B,B-s标量s分别与B元素之差B.A（一定与上相同） s.*A标量s分别与B元素之积exp(A)以自然数e为底，分别以A的元素为 指数，求幂 s./B,B.ss分别被B的元素除log(A)对A的各元素求对数 A.nA的每个元素自乘n次sqrt(A)对A的各元素求平方根 A.p对A各元素分别求非整数幂f(A)求各元素的函数值 p.A以p为底，分别以A的元素为 指数求幂 ABA、B阵对应元素间的关系运算 A+B对应元素相加ABA、B阵对应元素间的逻辑运算 运算指定形式 * 18 七、标准矩阵生成函数和数组操作函数 1通用的标准矩阵 常用的产生通用特殊矩阵的函数有： zeros：产生全0矩阵(零矩阵)。 ones：产生全1矩阵(幺矩阵)。 eye：产生单位矩阵。 rand：产生01间均匀分布的随机矩阵。 randn：产生均值为0，方差为1的标准正态分布随机矩阵。 * 19 【例】分别建立33、32和与矩阵A同样 大小的零矩阵。 解： (1) zeros(3) (2) zeros(3,2) (3) 设A为23矩阵，则可以用zeros(size(A) 建立一个与矩阵A同样大小零矩阵。 A=1 2 3;4 5 6; %产生一个23阶矩阵A zeros(size(A) %产生一个与矩阵A同样 大小的零矩阵 * 20 【例】建立随机矩阵： (1) 在区间20,50内均匀分布的5阶随机矩 阵。 (2) 均值为0.6、方差为0.1的5阶正态分布随 机矩阵。 解： x=20+(50-20)*rand(5) y=0.6+sqrt(0.1)*randn(5) * 21 2用于专门学科的特殊矩阵 (1)魔方矩阵 魔方矩阵有一个有趣的性质，其每行、每列及两条 对角线上的元素和都相等。对于n阶魔方阵，其元素 由1,2,3,n2共n2个整数组成。 MATLAB提供了求魔方矩阵的函数magic(n)，其 功能是生成一个n阶魔方阵。 （2）范得蒙矩阵vander(V) （3）希尔伯特矩阵hilb(n) 希尔伯特矩阵的逆矩阵invhilb(n) （4）托普利兹矩阵toeplitz(x) （5）伴随矩阵compan(p) （6）帕斯卡矩阵pascal(n)* 22 【例】将101125等25个数填入一个5行5列的表格中 ，使其每行每列及对角线的和均为565。 解：M=100+magic(5) M = 117 124 101 108 115 123 105 107 114 116 104 106 113 120 122 110 112 119 121 103 111 118 125 102 109 * 23 【例】标准数组产生的演示。 ones(1,2) %产生长度为2的全1行数组 ans = 1 1 ones(2) %产生（22）全1阵 ans = 1 1 1 1 randn(state,0) %把正态随机数发生器置0 randn(2,3) %产生（23）的正态随机阵 ans = -0.4326 0.1253 -1.1465 -1.6656 0.2877 1.1909 D=eye(3) %产生（33）的单位阵 D = 1 0 0 0 1 0 0 0 1 * 24 diag(D) %取D阵的对角元素 ans = 1 1 1 diag(diag(D) %内diag取D的对角元素，外diag利用一维 数组生成对角阵 ans = 1 0 0 0 1 0 0 0 1 repmat(D,1,3) %在水平方向“铺放”3个D阵 ans = Columns 1 through 9 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 * 25 3、数组操作函数 【例 】diag与reshape的使用演示。 a=-4:4 %产生一维数组 A=reshape(a,3,3) %把一维数组a重排成（33）的二维数组 a = Columns 1 through 9 -4 -3 -2 -1 0 1 2 3 4 A = -4 -1 2 -3 0 3 -2 1 4 a1=diag(A,1) %取A阵“第一上对角线”元素 a1 = -1 3 * 26 A1=diag(a1,-1) %产生以a1数组元素为“第一下对角线 ”元素的二维数组 A1 = 0 0 0 -1 0 0 0 3 0 A. %转置 ans = -4 -3 -2 -1 0 1 2 3 4 flipud(A) %上下对称交换 ans = -2 1 4 -3 0 3 -4 -1 2 * 27 八、数组构成技法综合 a=1 2; b=3 4; c=5;6; d=a;b; e=d c; f=e e;a b a; * 28 【例】数组的扩展。 （1）数组的赋值扩展法 A=reshape(1:9,3,3) %创建（33）数组A A = 1 4 7 2 5 8 3 6 9 A(5,5)=111 %扩展为（55）数组，扩展部分除(5,5)元素 为111外，其余均为0 A = 1 4 7 0 0 2 5 8 0 0 3 6 9 0 0 0 0 0 0 0 0 0 0 0 111 * 29 A(:,6)=222 %标量对子数组赋值，并扩展为（56）数组 A = 1 4 7 0 0 222 2 5 8 0 0 222 3 6 9 0 0 222 0 0 0 0 0 222 0 0 0 0 111 222 （2）多次寻访扩展法 AA=A(:,1:6,1:6) %相当于指令repmat(A,1,2) AA = 1 4 7 0 0 222 1 4 7 0 0 222 2 5 8 0 0 222 2 5 8 0 0 222 3 6 9 0 0 222 3 6 9 0 0 222 0 0 0 0 0 222 0 0 0 0 0 222 0 0 0 0 111 222 0 0 0 0 111 222 * 30 （3）合成扩展法 B=ones(2,6) %创建（26）全1数组 B = 1 1 1 1 1 1 1 1 1 1 1 1 AB_r=A;B %行数扩展合成 AB_r = 1 4 7 0 0 222 2 5 8 0 0 222 3 6 9 0 0 222 0 0 0 0 0 222 0 0 0 0 111 222 1 1 1 1 1 1 1 1 1 1 1 1 * 31 AB_c=A,B(:,1:5) %列数扩展合成 AB_c = 1 4 7 0 0 222 1 1 2 5 8 0 0 222 1 1 3 6 9 0 0 222 1 1 0 0 0 0 0 222 1 1 0 0 0 0 111 222 1 1 * 32 九、多项式的表达方式及其操作 1、多项式的创建 （1）多项式系数向量的直接输入法 （2）利用指令：P=poly(AR)产生多项式系数向量 若AR是方阵，则多项式P就是该方阵的特征多项式 若AR=ar1 ar2 ar3arn，则AR的元素被认为是 多项式P的根 【例 】求3阶方阵A的特征多项式。 A=11 12 13;14 15 16;17 18 19; PA=poly(A) PPA=poly2str(PA,s) PA = 1.0000 -45.0000 -18.0000 0.0000 PPA = s3 - 45 s2 - 18 s + 1.8303e-014 * 33 2、多项式运算函数 指 令含 义 p=conv(p1,p2)P是多项式p1和p2的乘积多项式 q,r=deconv(p1,p2)多项式p1被p2除的商多项式为q，而余多项式为r p=ploy(AR)求方阵AR的特征多项式p；或求向量AR指定根所对应的多 项式 dp=polyder(p)求多项式p的导数多项式dp dp=polyder(p1,p2)求多项式p1,p2乘积的导数多项式dp Num,Den=ployder(p1,p2)对有理分式(p1/p2)求导数所得的有理分式为(Num/Den) P=polyfit(x,y,n)求x,y向量给定的数据的n阶拟合多项式p pA=polyval(p,S)按数组运算规则计算多项式值；p为多项式，S为矩阵 PM=polyvalm(p,S)按矩阵运算规则计算多项式值；p为多项式，S为矩阵 r,p,k=residue(b,a)部分分式展开，b,a分别是分子、分母多项式系数向量；r 、p、k分别是留数、极点和直项 r=roots(p)求多项式p的根 * 34 【例】求 的“商”及“余”多项式。 p1=conv(1,0,2,conv(1,4,1,1); %计算分子多项式 p2=1 0 1 1; %注意缺项补零 q,r=deconv(p1,p2); cq=商多项式为 ; cr=余多项式为 ; disp(cq,poly2str(q,s),disp(cr,poly2str(r,s) 商多项式为 s + 5 余多项式为 5 s2 + 4 s + 3 * 35 十、多项式运算 1、多项式的创建 （1）多项式系数向量的直接输入法 （2）利用指令：P=poly(AR)产生多项式系数向量 若AR是方阵，则多项式P就是该方阵的特征多项式 若AR=ar1 ar2 ar3arn，则AR的元素被认为是多项式P的根 【例 】求3阶方阵A的特征多项式。 A=11 12 13;14 15 16;17 18 19; PA=poly(A) PPA=poly2str(PA,s) PA = 1.0000 -45.0000 -18.0000 0.0000 PPA = s3 - 45 s2 - 18 s + 1.8303e-014 * 36 2、多项式运算函数 指 令含 义 p=conv(p1,p2)P是多项式p1和p2的乘积多项式 q,r=deconv(p1,p2)多项式p1被p2除的商多项式为q，而余多项式为r p=ploy(AR)求方阵AR的特征多项式p；或求向量AR指定根所对应的多 项式 dp=polyder(p)求多项式p的导数多项式dp dp=polyder(p1,p2)求多项式p1,p2乘积的导数多项式dp Num,Den=ployder(p1,p2)对有理分式(p1/p2)求导数所得的有理分式为(Num/Den) P=polyfit(x,y,n)求x,y向量给定的数据的n阶拟合多项式p pA=polyval(p,S)按数组运算规则计算多项式值；p为多项式，S为矩阵 PM=polyvalm(p,S)按矩阵运算规则计算多项式值；p为多项式，S为矩阵 r,p,k=residue(b,a)部分分式展开，b,a分别是分子、分母多项式系数向量；r 、p、k分别是留数、极点和直项 r=roots(p)求多项式p的根 * 37 【例】求 的“商”及“余”多项式。 p1=conv(1,0,2,conv(1,4,1,1); %计算分子多项式 p2=1 0 1 1; %注意缺项补零 q,r=deconv(p1,p2); cq=商多项式为 ; cr=余多项式为 ; disp(cq,poly2str(q,s),disp(cr,poly2str(r,s) 商多项式为 s + 5 余多项式为 5 s2 + 4 s + 3 【例】 求多项式的微分polyder(p) polyder(a,b): 求多项式a,b乘积的微分 p,q=polyder(a,b): 求多项式a,b商的微分 a=1 2 3 4 5; poly2str(a,x) ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyder(a) b = 4 6 6 4 poly2str(b,x) ans =4 x3 + 6 x2 + 6 x + 4 * 38 3代数多项式求值 polyval函数用来求代数多项式的值，其调用格式为： Y=polyval(P,x) 若x为一数值，则求多项式在该点的值；若x为向量或 矩阵，则对向量或矩阵中的每个元素求其多项式的值 。 【例】已知多项式x4+8x3-10，分别取x=1.2和一个23矩 阵为自变量计算该多项式的值。 A=1,8,0,0,-10; %4次多项式系数 x=1.2; %取自变量为一数值 y1=polyval(A,x) y1 =5.8976 x=-1,1.2,-1.4;2,-1.8,1.6; %给出一个矩阵x y2=polyval(A,x) y2 =-17.0000 5.8976 -28.1104 70.0000 -46.1584 29.3216 * 39 4矩阵多项式求值 polyvalm函数用来求矩阵多项式的值，其调用格式与 polyval相同，但含义不同。polyvalm函数要求x为方阵, 它以方阵为自变量求多项式的值。设A为方阵，P代表 多项式x3-5x2+8，那么polyvalm(P,A)的含义是： A*A*A-5*A*A+8*eye(size(A) 而polyval(P,A)的含义是： A.*A.*A-5*A.*A+8*ones(size(A) 【例】仍以多项式x4+8x3-10为例，取一个22矩阵为自 变量分别用polyval和polyvalm计算该多项式的值。 A=1,8,0,0,-10; %多项式系数 x=-1,1.2;2,-1.8; %给出一个矩阵x * 40 y1=polyval(A,x) %计算代数多项式的值 y1 = -17.0000 5.8976 70.0000 -46.1584 y2=polyvalm(A,x) %计算矩阵多项式的值 y2 = -60.5840 50.6496 84.4160 -94.3504 * 41 5 多项式求根 n次多项式具有n个根，当然这些根可能是实根， 也可能含有若干对共轭复根。MATLAB提供的 roots函数用于求多项式的全部根，其调用格式为 ： x=roots(P) 其中P为多项式的系数向量，求得的根赋给向量x ，即x(1),x(2),x(n)分别代表多项式的n个根。 【例】求多项式x4+8x3-10的根。 命令如下： A=1,8,0,0,-10; x=roots(A) x = -8.0194 1.0344 -0.5075 + 0.9736i -0.5075 - 0.9736i * 42 十一、数据统计处理 (一)最大值和最小值 MATLAB提供的求数据序列的最大值和最小值的 函数分别为max和min，两个函数的调用格式和操 作过程类似。 1求向量的最大值和最小值 求一个向量X的最大值的函数有两种调用格式，分 别是： (1) y=max(X)：返回向量X的最大值存入y，如果X中包 含复数元素，则按模取最大值。 (2) y,I=max(X)：返回向量X的最大值存入y，最大值 的序号存入I，如果X中包含复数元素，则按模取 最大值。 * 43 2求矩阵的最大值和最小值 求矩阵A的最大值的函数有3种调用格式，分别是： (1) max(A)：返回一个行向量，向量的第i个元素是矩 阵A的第i列上的最大值。 (2) Y,U=max(A)：返回行向量Y和U，Y向量记录A 的每列的最大值，U向量记录每列最大值的行号。 (3) max(A,dim)：dim取1或2。dim取1时，该函数 和max(A)完全相同；dim取2时，该函数返回一个 列向量，其第i个元素是A矩阵的第i行上的最大值 。 * 44 3两个向量或矩阵对应元素的比较 函数max和min还能对两个同型的向量或矩阵进行 比较，调用格式为： (1) U=max(A,B)：A,B是两个同型的向量或矩阵，结果 U是与A,B同型的向量或矩阵，U的每个元素等于 A,B对应元素的较大者。 (2) U=max(A,n)：n是一个标量，结果U是与A同型的 向量或矩阵，U的每个元素等于A对应元素和n中的 较大者。 min函数的用法和max完全相同。 * 45 (二)求和与求积 数据序列求和与求积的函数是sum和prod，其使用方法 类似。设X是一个向量，A是一个矩阵，函数的调用格 式为： sum(X)：返回向量X各元素的和。 prod(X)：返回向量X各元素的乘积。 sum(A)：返回一个行向量，第i个元素是A的第i列的元素和 。 prod(A)：返回一个行向量，第i个元素是A的第i列的元素乘 积。 sum(A,dim)：dim=1时，该函数等同于sum(A)；dim=2时， 返回一个列向量，其第i个元素是A的第i行的各元素之 和。 prod(A,dim)：dim=1时，该函数等同于prod(A)；dim=2时 ，返回一个列向量，其第i个元素是A的第i行的各元素 乘积。 * 46 (三)平均值和中值 求数据序列平均值的函数是mean，求数据序列中值 的函数是median。两个函数的调用格式为： mean(X)：返回向量X的算术平均值。 median(X)：返回向量X的中值。 mean(A)：返回一个行向量，其第i个元素是A的第i列 的算术平均值。 median(A)：返回一个行向量，其第i个元素是A的第i列 的中值。 mean(A,dim)：dim=1时，该函数等同于mean(A)； dim=2时，返回一个列向量，其第i个元素是A的第i行 的算术平均值。 median(A,dim)：dim=1时，该函数等同于median(A)； dim=2时，返回一个列向量，其第i个元素是A的第i行 的中值。 * 47 (四)累加和与累乘积 在MATLAB中，使用cumsum和cumprod函数能方便地求得 向量和矩阵元素的累加和与累乘积向量，函数的调用格式为 ： cumsum(X)：返回向量X累加和向量。 cumprod(X)：返回向量X累乘积向量。 cumsum(A)：返回一个矩阵，其第i列是A的第i列的累加和向 量。 cumprod(A)：返回一个矩阵，其第i列是A的第i列的累乘积向 量。 cumsum(A,dim)：dim=1时，函数等同于cumsum(A)；dim=2 时，返回一个矩阵，其第i行是A的第i行的累加和向量。 cumprod(A,dim)：dim=1时，函数等同于cumprod(A)；dim=2 时，返回一个向量，其第i行是A的第i行的累乘积向量。 * 48 (五) 标准方差与相关系数 1求标准方差 在MATLAB中，提供了计算数据序列的标准方差的函数 std。 std(X) ：返回一个标准方差。 std(A)：返回一个行向量，每个元素是矩阵A各列或各行 的标准方差。 Y=std(A,flag,dim) 缺省flag=0，dim=1 当dim=1时，求各列元素的标准方差； 当dim=2时，则求各行元素的标准方差。 当flag=0时，按S1所列公式计算标准方差； 当flag=1时，按S2所列公式计算标准方差。 * 49 【例】对二维矩阵x，从不同维方向求出其标准方差。 命令如下： x=4,5,6;1,4,8; y1=std(x,0,1) y2=std(x,1,1) y3=std(x,0,2) y4=std(x,1,2) 2相关系数 MATLAB提供了corrcoef函数，可以求出数据的相关系 数矩阵。corrcoef函数的调用格式为： corrcoef(X)：返回从矩阵X形成的一个相关系数矩阵。此 相关系数矩阵的大小与矩阵X一样。它把矩阵X的每列 作为一个变量，然后求它们的相关系数。 corrcoef(X,Y)：在这里，X,Y是向量，它们与 corrcoef(X,Y)的作用一样。 * 50 【例】生成满足正态分布的100005随机矩阵，然后 求各列元素的均值和标准方差，再求这5列随机数 据的相关系数矩阵。 命令如下： X=randn(10000,5); M=mean(X) D=std(X) R=corrcoef(X) * 51 (六)排序 MATLAB中对向量X是排序函数是sort(X)，函数返回 一个对X中的元素按升序排列的新向量。 Y,I=sort(A,dim)：对矩阵A的各列或各行重新排序 若dim=1，则按列排； 若dim=2，则按行排。 Y是排序后的矩阵，而I记录Y中的元素在A中位置。 【例】对下列矩阵做各种排序。 A=1,-8,5;4,12,6;13,7,-13; sort(A) %对A的每列按升序排序 -sort(-A,2) %对A的每行按降序排序 X,I=sort(A) %对A按列排序，并将每个元素所在的 行号送矩阵I * 52 (一) 一维数据插值 数据插值的任务是根据采集到的离散数据构造一个函 数g(x)，既与真实函数f(x)接近，又有很好的性质。 Y1=interp1(X,Y,X1,method) 一维插值函数根据X,Y的值，计算函数在X1处的值。 X,Y是两个等长的已知向量，分别描述采样点和样本值 ，X1是一个向量或标量，描述欲插值的点，Y1是与X1 等长的插值结果。method是插值方法，例如： linear线性插值 nearest最近点插值 cubic三次多项式插值 spline三次样条插值 或Y1=spline(X,Y,X1)计算三次样条插值的专门函数 十二、数据插值 * 53 注意：X1的取值范围不能超出X的给定范围，否则，会 给出“NaN”错误。 【例】某观测站测得某日6:00时至18:00时之间每隔2小时 的室内外温度()，用3次样条插值分别求得该日室内外 6:30至17:30时之间每隔2小时各点的近似温度()。 解：设时间变量h为一行向量，温度变量t为一个两列矩阵 ，其中第一列存放室内温度，第二列储存室外温度。命 令如下： h =6:2:18; t=18,20,22,25,30,28,24;15,19,24,28,34,32,30; XI =6.5:2:17.5 YI=interp1(h,t,XI, spline) %用3次样条插值计算 * 54 (二)二维数据插值 在MATLAB中，提供了解决二维插值问题的函数 interp2，其调用格式为： Z1=interp2(X,Y,Z,X1,Y1,method) 其中X,Y是两个向量，分别描述两个参数的采样点，Z 是与参数采样点对应的函数值，X1,Y1是两个向量或 标量，描述欲插值的点。Z1是根据相应的插值方法得 到的插值结果。 method的取值与一维插值函数相同。 X,Y,Z也可以是矩阵形式。 同样，X1,Y1的取值范围不能超出X,Y的给定范围，否 则，会给出“NaN”错误。 * 55 【例】设z=x2+y2，对z函数在0,10,2区域内进行插值 。 x=0:0.1:1; y=0:0.2:2; X,Y=meshgrid(x,y); %产生自变量网格坐标 Z=X.2+Y.2; %求对应的函数值 interp2(x,y,Z,0.5,0.5) %在(0.5,0.5)点插值 【例】某实验对一根长10米的钢轨进行热源的温度传播 测试。用x表示测量点0:2.5:10(米)，用h表示测量时 间0:30:60(秒)，用T表示测试所得各点的温度()。 试用线性插值求出在一分钟内每隔20秒、钢轨每隔1 米处的温度TI。 x=0:2.5:10;h=0:30:60; T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41; xi=0:10;hi=0:20:60; TI=interp2(x,h,T,xi,hi);mesh(xi,hi,TI) * 56 在MATLAB中，用polyfit函数来求得最小二乘拟合多项 式的系数，再用polyval函数按所得的多项式计算所给出 的点上的函数近似值。 P,S=polyfit(X,Y,m) 函数根据采样点X和采样点函数值Y，产生一个m次多项 式P及其在采样点的误差向量S。其中X,Y是两个等长的 向量，P是一个长度为m+1的向量，P的元素为多项式系 数。 polyval函数的功能是按多项式的系数计算x点多项式的值 。 十三、曲线拟合 【例】用一个三次多项式在区间0,2内逼近函数sin(x)。 X=linspace(0,2*pi,50); Y=sin(X); P,S=polyfit(X,Y,3) %得到3次多项式的系数和误差 * 57 * 58/113 function spfit x=0 1 2 3 4 5; y=5 2 1 1 2 3; fx=polyfit(x,y,2) %对数据进行二次拟合,输出结果 yvalue=polyval(fx,x); ymean=mean(y); Rsquare=(norm(yvalue-ymean)/norm(y- ymean)2 %求R的平方并输出 x1=0:0.2:5; y1=polyval(fx,x1); plot(x,y,b*,x1,y1,g) * 59/113 十四、数值积分与微分 (一)数值积分 1、变步长辛普生法 基本思想都是将整个积分区间a,b分成n个子区间 xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分 问题就分解为求和问题。 MATLAB给出了quad和quadl函数来求定积分。 I,n=quad(filename,a,b,tol,trace) I,n=quadl(filename,a,b,tol,trace) 其中filename是被积函数名。a和b分别是定积分的下限 和上限。tol用来控制积分精度，缺省时取tol=10-6。 trace控制是否展现积分过程，若取非0则展现积分过程 ，取0则不展现，缺省时取trace=0。返回参数I即定积分 值，n为被积函数的调用次数。 * 60 【例】求定积分 解：(1) 建立被积函数文件fesin.m： function f=fesin(x) f=exp(-0.5*x).*sin(x+pi/6); (2) 调用数值积分函数quad求定积分： S,n=quad(fesin,0,3*pi) 或S,n=quadl(fesin,0,3*pi) 或S,n=quad(fesin,0,3*pi) S = 0.9008 n = 77 quad函数采用自适应变步长方法求取定积分 quadl函数则采用Lobbato方法求取定积分，计算精 度和速度要高于quad函数方法。利用quadgk函数求 解广义积分问题。 * 61 2被积函数由一个表格定义 在MATLAB中，对由表格形式定义的函数关系的求定 积分问题用trapz(X,Y)函数。其中向量X,Y定义函数关 系Y=f(X)。 【例】用trapz函数计算定积分 命令如下： X=1:0.01:2.5; Y=exp(-X); %生成函数关系数据向量 trapz(X,Y) ans = 0.2858 * 62 3、二重积分的数值求解 使用MATLAB提供的dblquad函数就可以直接求出 二重定积分的数值解。该函数的调用格式为： I=dblquad(f(x,y),a,b,c,d,tol,trace) 该函数求f(x,y)在a,bc,d区域上的二重定积分。 参数tol，trace的用法与函数quad完全相同。 【例】计算二重定积分 解：(1) 建立一个函数文件fxy.m： function f=fxy(x,y) global ki; ki=ki+1; %ki用于统计被积函数的调用次数 f=exp(-x.2/2).*sin(x.2+y); (2)调用dblquad函数 global ki;ki=0; I=dblquad(fxy,-2,2,-1,1) I = 1.5745 * 63 4、一般二重积分的数值求解 MATLAB中没有提供求解一般的双重积分问题的函数 在Mathworks公司的网站上有一个数值积分工具箱NIT （Numerical Integration Toolbox）可以解决此问题。 利用工具箱中的函数quad2dggen()可以直接求解此类 问题。 J=quad2dggen(Fun,Flower,Fupper,xm,xM,tol) 默认tol=10-3，Fun描述二元被积函数，需要注意是： 该函数指定的积分顺序是先x后y，若要求先y后x的积分 ，需要交换被积函数中x和y的顺序。 * 64 【例16】计算二重定积分 解： fh=(x)sqrt(1-x.2/2); fl=(x)-sqrt(1-x.2/2); f=(y,x)exp(-x.2/2).*sin(x.2+y); y=quad2dggen(f,fl,fh,-1/2,1,eps) y = 0.4119 * 65 5、三重及多重积分的数值求解 长方体区域的三重积分问题的函数 采用triplequad()函数可以求解此类问题。 J=triplequad(Fun,xm,xM,ym,yM,zm,zM,tol,quadl) 默认tol=10-6。 NIT工具箱也可以解决多重超长方体边界的定积分问题 。使用quadndg()函数。 该问题的求解语句为： J=quadndg(Funx1m,x2m,xpm,x1M,x2M,xpM,tol)* 66 【例】计算三重定积分 解： f=(x,y,z)4*x.*z.*exp(-x.*x.*y-z.*z); triplequad(f,0,2,0,pi,0,pi,1e-7,equadl) ans = 3.1081 该工具箱中还有单重积分函数quadg()的调用格式同 quad()，其效率也高于quadl()。 * 67 (二)数值微分 1、数值差分与差商 高等数学关心的是导函数的形式和性质，而数值分 析关心的问题是怎样计算导函数： f (x)=g(x) 在一串离散点X=(x1,x2,xn)的近似值 G=(g1,g2,gn)以及所计算的近似值有多大误差。 引进记号： 称 分别为函数在x点处以 h(h0)为步长的向前差分、向后差分和中心差分。 * 68 称 分别为函数在x点处 以h(h0)为步长的向前差商、向后差商和中心差商。 当步长h(h0)充分小时，函数f在点x的微分接近于函数 在该点的任意种差分。 f在点x的导数接近于函数在该点 的任意种差商。 2、数值微分的实现 在MATLAB中，没有直接提供求数值导数的函数， 只有计算向前差分的函数diff，其调用格式为： DX=diff(X) ：计算X的向前差分