第三章数值数组及向量化运算(matlab基础教程)

第三章数值数组及向量化运算(matlab基础教程)

第3章 数值数组及向量化运算

数值数组（Numeric Array）和数组运算（Array Operations）始终是MATLAB的核心内容。本书从第3章起，全部注意力将集中于数值数组及其运算。 本章系统阐述：数组浮点算法的特点；一、二维数值数组的创建和寻访；数组运算和向量化编程；实现数组运算的基本函数；常用标准数组生成函数和数组构作技法；非数NaN、“空”数组概念和应用；关系和逻辑操作。

3.1 数值计算的特点和地位

【例3.1-1】已知f(t) t2cost，求s(x) （1）符号计算解法

syms t x ft=t^2*cos(t) sx=int(ft,t,0,x) ft =

t^2*cos(t) sx =

x^2*sin(x)-2*sin(x)+2*x*cos(x)

x

f(t) dt。

（2）数值计算解法

dt=0.05; t=0:dt:5; Ft=t.^2.*cos(t); Sx=dt*cumtrapz(Ft); t(end-4:end)

Sx(end-4:end)

plot(t,Sx,'.k','MarkerSize',12) xlabel('x'),ylabel('Sx'),grid on

ans =

4.8000 4.8500 4.9000 4.9500 5.0000

ans =

-20.1144 -19.9833 -19.7907 -19.5345 -19.2131

图 3.1-1 在区间[0, 5]采样点上算得的定积分值

【例3.1-2】已知

f(t) e sin(t)，求s(x) 0f(t) dt。

4

（1）符号计算解法

syms t x

ft=exp(-sin(t)) sx=int(ft,t,0,4) ft =

exp(-sin(t))

Warning: Explicit integral could not be found. In at 58 sx =

int(exp(-sin(t)),t = 0 .. 4)

（2）数值计算解法

dt=0.05; t=0:dt:4; Ft=exp(-sin(t)); Sx=dt*cumtrapz(Ft); Sx(end)

plot(t,Ft,'*r','MarkerSize',4) hold on

plot(t,Sx,'.k','MarkerSize',15) hold off xlabel('x')

legend('Ft','Sx') ans =

3.0632

图 3.1-2 在区间[0, 4]中间的被积函数及其原函数的离散计算结果

3.2

3.2.1 一 二

数值数组的创建和寻访

一维数组的创建 递增/减型一维数组的创建 通用型一维数组的创建

【例3.2-1】一维数组的常用创建方法举例。

a1=1:6 a2=0:pi/4:pi a3=1:-0.1:0 a1 =

1 2 3 4 5 6 a2 =

0 0.7854 1.5708 2.3562 3.1416

a3 =

Columns 1 through 6

1.0000 0.9000 0.8000 0.7000 0.6000 0.5000 Columns 7 through 11

0.4000 0.3000 0.20XX年 0.1000 0

b1=linspace(0,pi,4) b2=logspace(0,3,4) b1 =

0 1.0472 2.0944 3.1416 b2 =

1 10 100 1000

c1=[2 pi/2 sqrt(3) 3+5i] c1 =

Columns 1 through 3

2.0000 1.5708 1.7321 Column 4

3.0000 + 5.0000i

rand('state',0) c2=rand(1,5) c2 =

0.9501 0.2311 0.6068 0.4860 0.8913

3.2.2 一

二维数组的创建 小规模数组的直接输入法

【例3.2-2】在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

二 中规模数组的数组编辑器创建法

【例3.2-3】根据现有数据创建一个(3 8)的数组。

图3.2-1 利用数组编辑器创建中规模数组

三 中规模数组的M文件创建法

【例3.2-4】创建和保存数组 AM的 MyMatrix.m 文件。

（1） 打开文件编辑调试器，并在空白填写框中输入所需数组（见图3.2-2）。

（2） 最好，在文件的首行，编写文件名和简短说明，以便查阅（见图3.2-2）。 （3） 保存此文件，并且文件起名为 MyMatrix.m 。

（4） 以后只要在MATLAB指令窗中，运行 MyMatrix.m文件 ，数组 AM 就会自动生成

于 MATLAB内存中。

图3.2-2 利用M文件创建数组

四 利用MATLAB函数创建数组

【例3.2-5】标准数组产生的演示。

ones(2,4) ans =

1 1 1 1 1 1 1 1

randn('state',0) randn(2,3) ans =

-0.4326 0.1253 -1.1465 -1.6656 0.2877 1.1909

D=eye(3) D =

1 0 0 0 1 0 0 0 1

diag(D) ans = 1 1 1

diag(diag(D)) ans =

1 0 0 0 1 0 0 0 1

randsrc(3,20,[-3,-1,1,3],1) ans =

Columns 1 through 13

3 1 3 -1 -3 3 -3 -3 -1 3 -1 -1 -3 1 3 1 3 -1 1 1 1 1 1 1 1 3 3 -1 -3 -1

1 -1 3 -1 -1 1 1 -1 -3 Columns 14 through 20

1 1 -3 3 -1 1 3 -1 -1 -3 -1 -1 1 -3 3 -1 -1 3 -1 3 3

3.2.3 二维数组元素的标识和寻访

【例3.2-6】本例演示：数组元素及子数组的各种标识和寻访格式；冒号的使用；end的作用。

A=zeros(2,6) A(:)=1:12 A =

0 0 0 0 0 0 0 0 0 0 0 0 A =

1 3 5 7 9 11 2 4 6 8 10 12

A(2,4) A(8) ans = 8 ans = 8

A(:,[1,3]) A([1,2,5,6]') ans =

1 5 2 6 ans = 1 2 5 6

A(:,4:end) ans =

7 9 11 8 10 12

A(2,1:2:5)=[-1,-3,-5] A =

1 3 5 7 9 11 -1 4 -3 8 -5 12

B=A([1,2,2,2],[1,3,5] ) B =

1 5 9 -1 -3 -5 -1 -3 -5 -1 -3 -5

L=A3 A(L)=NaN L =

1 0 0 0 0 0 1 0 1 0 1 0 A =

NaN 3 5 7 9 11 NaN 4 NaN 8 NaN 12

3.2.4

a=1:8

数组构作技法综合

【例 3.2-7】数组操作函数reshape, diag, repmat的用法；空阵 [ ]

删除子数组的用法。

A=reshape(a,4,2) A=reshape(A,2,4) a =

1 2 3 4 5 6 7 8 A =

1 5 2 6 3 7 4 8 A =

1 3 5 7 2 4 6 8

b=diag(A) B=diag(b) b = 1 4 B =

1 0 0 4

D1=repmat(B,2,4) D1 =

1 0 1 0 1 0 1 0 0 4 0 4 0 4 0 4 1 0 1 0 1 0 1 0 0 4 0 4 0 4 0 4

D1([1,3],: )=[ ] D1 =

0 4 0 4 0 4 0 4 0 4 0 4 0 4 0 4

【例3.2-8】函数flipud, fliplr, rot90对数组的操作体现着“矩阵变

换”。

A=reshape(1:9,3,3) A =

1 4 7 2 5 8 3 6 9

B=flipud(A) B =

3 6 9 2 5 8 1 4 7

C=fliplr(A) C =

7 4 1 8 5 2 9 6 3

D=rot90(A,2) D =

9 6 3 8 5 2 7 4 1

3.3 数组运算

3.3.1 一 二 三 四 3.3.2

数组运算的由来和规则 函数关系数值计算模型的分类 提高程序执行性能的三大措施 数组运算规则

数组运算符及数组运算函数

数组运算和向量化编程

u

，其中r, u, i分别是电阻（欧姆）、电压（伏特）、电流（安i

【例 3.3-1】欧姆定律：r

培）。验证实验：据电阻两端施加的电压，测量电阻中流过的电流，然后据测得的电压、电流计算平均电阻值。（测得的电压电流具体数据见下列程序）。 （1）非向量化程序

clear

vr=[0.89, 1.20, 3.09, 4.27, 3.62, 7.71, 8.99, 7.92, 9.70, 10.41];

ir=[0.028, 0.040, 0.100, 0.145, 0.118, 0.258, 0.299, 0.257, 0.308,

0.345]; % -------------------- L=length(vr); for k=1:L r(k)=vr(k)/ir(k); end

% --------------------------- sr=0;

for k=1:L sr=sr+r(k); end

rm=sr/L rm =

30.5247

（2）向量化程序

clear

vr=[0.89, 1.20, 3.09, 4.27, 3.62, 7.71, 8.99, 7.92, 9.70, 10.41];

ir=[0.028, 0.040, 0.100, 0.145, 0.118, 0.258, 0.299, 0.257, 0.308, 0.345];

r=vr./ir rm=mean(r) r =

Columns 1 through 7

31.7857 30.0000 30.9000 29.4483 30.6780 29.8837 30.0669 Columns

8 through 10

30.8171 31.4935 30.1739 rm =

30.5247

【例 3.3-2】用间距为0.1的水平线和垂直线均匀分割x [ 5,5],y

[ 2.5,2.5]的矩形域，在所有水平线和垂直线交点上计算函数z sin|xy|的值，并图示。 （1）

clear

x=-5:0.1:5;

y=(-2.5:0.1:2.5)';

N=length(x); M=length(y); for ii=1:M

for jj=1:N

X0(ii,jj)=x(jj); Y0(ii,jj)=y(ii);

Z0(ii,jj)=sin(abs(x(jj)*y(ii))); end end

（2）

[X,Y]=meshgrid(x,y); Z=sin(abs(X.*Y));

（3）

norm(Z-Z0) ans = 0

（4）

surf(X,Y,Z) xlabel('x') ylabel('y') shading interp

view([190,70])

图 3.3-1 指定域上的二元函数图形

3.4

3.4.1

“非数”和“空”数组

非数NaN

【例3.4-1】非数的产生和性质演示。 （1）非数的产生

a=0/0,b=0*log(0),c=inf-inf

a = NaN b = NaN c =

NaN

（2）非数的传递性

0*a,sin(a)

ans = NaN ans = NaN

（3）非数的属性判断

class(a) isnan(a) ans = double ans = 1

【例3.4-2】非数元素的寻访。

rand('state',0)

R=rand(2,5);R(1,5)=NaN;R(2,3)=NaN R =

0.9501 0.6068 0.8913 0.4565 NaN 0.2311 0.4860 NaN 0.0185 0.4447

LR=isnan(R) LR =

0 0 0 0 1 0 0 1 0 0

si=find(LR) [ri,ci]=ind2sub(size(R),si) [rj,cj]=find(LR)

disp('非数在二维数组R中的位置')

disp(['单下标时的第',int2str(si(1)),'和第',int2str(si(2)),'个元素']) si = 6 9

ri = 2 1 ci = 3 5 rj = 2 1 cj = 3 5

非数在二维数组R中的位置

单下标时的第6和第9个元素

3.4.2 “空”数组

【例3.4-3】关于“空”数组的算例。 （1）创建“空”数组的几种方法

a=[]

b=ones(2,0),c=zeros(2,0),d=eye(2,0) f=rand(2,3,0,4) a =

[] b =

Empty matrix: 2-by-0 c =

Empty matrix: 2-by-0 d =

Empty matrix: 2-by-0 f =

Empty array: 2-by-3-by-0-by-4

（2）“空”数组的属性

class(a) isnumeric(a) isempty(a) ans = double ans = 1 ans = 1

which a ndims(a) size(a) a is a variable. ans = 2 ans =

0 0

（3）“空”数组用于子数组的删除和大数组的大小收缩 A=reshape(-4:5,2,5)

A =

-4 -2 0 2 4 -3 -1 1 3 5

A(:,[2,4])=[] A =

-4 0 4 -3 1 5

3.5

关系操作和逻辑操作

3.5.1

关系操作

【例3.5-1】关系运算示例。

A=1:9,B=10-A

r0=(A4) r1=(A==B) A =

1 2 3 4 5 6 7 8 9 B =

9 8 7 6 5 4 3 2 1 r0 =

1 1 1 0 0 0 0 0 0 r1 =

0 0 0 0 1 0 0 0 0

【例3.5-2】关系运算应用。

t=-3*pi:pi/10:3*pi; y=sin(t)./t; tt=t+(t==0)*eps; yy=sin(tt)./tt;

subplot(1,2,1),plot(t,y),axis([-9,9,-0.5,1.2]), xlabel('t'),ylabel('y'),title('残缺图形')

subplot(1,2,2),plot(tt,yy),axis([-9,9,-0.5,1.2])

xlabel('tt'),ylabel('yy'),title('正确图形')

图3.5-1 采用近似极限处理前后的图形对照

3.5.2 逻辑操作

【例3.5-3】逻辑操作和关系操作。本例演示：逻辑、关系操作的组合；xor的作用。 （1）逻辑、关系操作的组合

A=[-2,-1,0,0,1,2,3]

L1=~(A1) %判断A中，哪些元素不大于1

L2=(A0)(A2) %判断A中，哪些元素大于0且小于3 A =

-2 -1 0 0 1 2 3 L1 =

1 1 1 1 1 0 0 L2 =

0 0 0 0 1 0 0

（2）xor的作用

A,B=[0,-1,1,0,1,-2,-3]

C=xor(A,B) %当A, B数组中，两个对应元素中仅一个为0时，给出1。否则为0。

【例3.5-4】试绘制如图3.5-2最下那幅子图所示的“正弦波sint的削顶半波整流波形”，削顶发生在每个周期的[60, 120]之间。

clear,t=linspace(0,3*pi,500);y=sin(t); z1=((tpi)|(t2*pi)).*y;

w=(tpi/3t2*pi/3)+(t7*pi/3t8*pi/3); wn=~w;

z2=w*sin(pi/3)+wn.*z1;

subplot(4,1,1),plot(t,y,':r'),axis([0,10,-1.5,1.5]) ylabel('y'),grid on

subplot(4,1,2),plot(t,z1,':r'),axis([0,10,-0.2,1.5]),ylabel('z1')

subplot(4,1,3),plot(t,wn,':r'),axis([0,10,-0.2,1.5]),ylabel('wn')

subplot(4,1,4),plot(t,z2,'-b'),axis([0,10,-0.2,1.5]),ylabel('z2') xlabel('t')

3.5.3 常用逻辑函数

习题3

1． 要求在闭区间[0,2 ]上产生具有10个等距采样点的一维数组。试用两种不同的指令实

现。

2． 由指令rand('state',0),A=rand(3,5)生成二维数组A，试求该数组中所有大于0.5的元素的

位置，分别求出它们的“全下标”和“单下标”。 〖答案〗

大于0.5的元素的全下标

行号 1 3 2 3 3 2 3 1 2 列号 1 1 2 2 3 4 4 5 5

大于0.5的元素的单下标

1 3 5 6 9 11 12 13 14

3． 在使用123作为rand随机数发生器的初始化状态的情况下，写

出产生长度为1000的“等

概率双位（即取-1，+1）取值的随机码”程序指令，并给出 -1码的数目。 〖答案〗 Na =

490

4． 已知矩阵A

12

，运行指令B1=A.^(0.5), B2=A^(0.5), 可以观察到不同运算方法所 34

得结果不同。（1）请分别写出根据B1, B2恢复原矩阵A的程序。（2）用指令检验所得的两个恢复矩阵是否相等。

5． 在时间区间 [0,10]中，绘制y 1 e 0.5tcos2t曲线。要求分别采取“标量循环运算法”

和“数组运算法”编写两段程序绘图。

6． 先运行clear,format long,rand('state',1),A=rand(3,3)，然后根据A写出两个矩阵：一个对

角阵B，其相应元素由A的对角元素构成；另一个矩阵C，其对角元素全为0，而其余元素与对应的A阵元素相同。

7． 先运行指令x=-3*pi:pi/15:3*pi; y=x; [X,Y]=meshgrid(x,y); warning

off;

Z=sin(X).*sin(Y)./X./Y; 产生矩阵Z。（1）请问矩阵Z中有多少个“非数”数据？（2）用指令surf(X,Y,Z); shading interp观察所绘的图形。（3）请写出绘制相应的“无裂缝”图形的全部指令。

〖答案〗

NumOfNaN = 181

1k 1 9k 1 2k 2 9k 2

，8． 下面有一段程序，企图用来解决如下计算任务：有矩阵Ak

k2k 10k

当k依次取10, 9, 8, 7, 6, 5, 4, 3, 2, 1时，计算矩阵Ak“各列元素的和”，并把此求和

14 28 结果存放为矩阵Sa的第k行。例如k 3时，A阵为25 29，此时它各列元 36 30

素 的和是一个(1 10)行数组 615 87 ，并把它保存为Sa的第3行。问题：该

段程序的计算结果对吗？假如计算结果不正确，请指出错误发生的根源，并改正之。 for k=10:-1:1

A=reshape(1:10*k,k,10); Sa(k,:)=sum(A); end Sa