ch3 数值阵列及向量化运算.docx
《ch3 数值阵列及向量化运算.docx》由会员分享，可在线阅读，更多相关《ch3 数值阵列及向量化运算.docx（31页珍藏版）》请在冰豆网上搜索。
ch3数值阵列及向量化运算
第3章数值阵列及向量化运算
.1数值计算的特点和地位
符号计算的局限性：
1）有很多问题无法解，2）求解时间过长
数值计算：
适用范围广，能处理各种复杂的函数关系，计算速度快，容量大。
【例3.1-1】已知
，求
。
（1）符号计算解法
symstx
ft=t^2*cos（t）
sx=int（ft,t,0,x）
ft=
t^2*cos（t）
sx=
sin（x）*（x^2-2）+2*x*cos（x）
（2）数值计算解法
dt=0.05;
t=0:
dt:
5;
Ft=t.^2.*cos（t）;
Sx=dt*cumtrapz（Ft）;%由一个个宽度为dt的小梯形面积的累加求Ft曲线下的面积
t（end-4:
end）
Sx（end-4:
end）
plot（t,Sx,'.k','MarkerSize',12）
xlabel（'x'）,ylabel（'Sx'）,gridon
ans=
4.80004.85004.90004.95005.0000
ans=
-20.1144-19.9833-19.7907-19.5345-19.2131
图3.1-1在区间[0,5]采样点上算得的定积分值
【例3.1-2】已知
，求
。
（1）符号计算解法
symstx
ft=exp（-sin（t））
sx=int（ft,t,0,4）
ft=
exp（-sin（t））
Warning:
Explicitintegralcouldnotbefound.
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）
holdon
plot（t,Sx,'.k','MarkerSize',15）
holdoff
xlabel（'x'）
legend（'Ft','Sx'）
补充：
cumtrapz是计算步长为dt的梯形的面积，得到一组和t一样长的面积数列，其算法可以是dt*cumptrapz（ft）,也可以是cumtrapz（t,ft）；若仅仅是只要输出最终积分结果，则可以用trapz（t,ft）。
ans=
3.0632
.2数值数组的创建和寻访
.2.1一维数组的创建
x=[1,3,5,7,9]逐个元素输入法
x=a:
inc:
b冒号生成法，inc缺省时步长为1
x=linspace（a,b,n）线性定点法
x=logspace（a,b,n）对数定点法，以底数为十的对数值等距取点，可以在缩短数值之间的差距，便于表示。
运用标准数组生成函数。
【例3.2-1】一维数组的常用创建方法举例。
a1=1:
6
a2=0:
pi/4:
pi
a3=1:
-0.1:
0
a1=
123456
a2=
00.78541.57082.35623.1416
a3=
Columns1through6
1.00000.90000.80000.70000.60000.5000
Columns7through11
0.40000.30000.20000.10000
b1=linspace（0,pi,4）
b2=logspace（0,3,4）%创建数组[100101102103]
b1=
01.04722.09443.1416
b2=
1101001000
c1=[2pi/2sqrt（3）3+5i]
c1=
Columns1through3
2.00001.57081.7321
Column4
3.0000+5.0000i
rand（'state',0）%利用标准数组生成函数产生均匀分布随机数组
c2=rand（1,5）
c2=
0.95010.23110.60680.48600.8913
补充：
rand（‘state’,i）表示以后具有同样声明的随机函数产生的随机函数一样，其中i用于标出第几个随机数组。
.2.2二维数组的创建
10一小规模数组的直接输入法
【例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.00005.4716+0.4177i0.6909
0.70714.82443.5000+1.0000i
三个要素：
整个输入数组必须以方括号“[]”为其首尾；
行与行间“;”或“Enter”；
同行中元素间“,”或“空格”。
10二中规模数组的数组编辑器创建法
【例3.2-3】根据现有数据创建一个
的数组。
图3.2-1利用数组编辑器创建中规模数组
10三中规模数组的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文件创建数组
补充：
注意分号的使用，若使用分号，则抑制结果的表达，只显示在工作空间中。
10四利用MATLAB函数创建数组
【例3.2-5】利用最常用标准数组生成函数产生标准数组的演示。
ones（2,4）%产生（2×4）全1数组
ans=
1111
1111
randn（'state',0）%把正态随机数发生器置0
randn（2,3）%产生正态随机阵
ans=
-0.43260.1253-1.1465
-1.66560.28771.1909
D=eye（3）%产生3×3的单位阵
D=
100
010
001
diag（D）%取D阵的对角元
ans=
1
1
1
diag（diag（D））%外diag利用一维数组生成对角阵
ans=
100
010
001
randsrc（3,20,[-3,-1,1,3],1）%在[-3,-1,1,3]上产生3×20均布随机数组，随机发生器的状态设置为1
ans=
Columns1through13
313-1-33-3-3-13-1-1-3
1313-111111113
3-1-3-11-13-1-111-1-3
Columns14through20
11-33-113
-1-1-3-1-11-3
3-1-13-133
.2.3二维数组元素的编址和寻访
【例3.2-6】本例演示：
数组元素及子数组的各种标识和寻访格式；冒号的使用；end的作用。
A=zeros（2,6）
A（:
）=1:
12%单下标法：
单下标全元素寻访
%补充：
此处还可以用A（:
）=1：
2：
23来构成一个等差数列。
A=
000000
000000
A=
1357911
24681012
A（2,4）%全下标法：
指定行、指定列
A（8）%单下标法：
单下标寻访
ans=
8
ans=
8
A（:
[1,3]）%全下标法：
全部行、指定列
%注意，此处使用的是中括号，在用subs函数时可以用花括号，也可以用中括号。
A（[1,2,5,6]'）%单下标法:
生成指定的一维行（或列）数组
ans=
15
26
ans=
1
2
5
6
补充：
单引号表示转置。
A（:
4:
end）%全下标法：
全部行、指定列，end表示最后一列。
ans=
7911
81012
A（2,1:
2:
5）=[-1,-3,-5]%全下标法：
指定行、指定列
A=
1357911
-14-38-512
B=A（[1,2,2,2],[1,3,5]）%全下标法：
指定行、指定列
B=
159
-1-3-5
-1-3-5
-1-3-5
补充：
行与列的结合为乘法，即每取一行，都要将后面对应的列取完。
相互转换指令：
ind2sub：
把单序号编址转换成全下标编址
sub2ind：
把全下标编址转换成到单序号编址
.2.4数组构作技法综合
【例3.2-7】数组操作函数reshape,diag,repmat的用法；空阵[]删除子数组的用法。
a=1:
8
A=reshape（a,4,2）
A=reshape（A,2,4）%改变行数和列数
a=
12345678
A=
15
26
37
48
A=
1357
2468
b=diag（A）%提取对角元素，。
B=diag（b）%生成对角阵
b=
1
4
B=
10
04
补充：
双重取对角元数即生成对角阵。
D1=repmat（B,2,4）%排列B模块
D1=
10101010
04040404
10101010
04040404
D1（[1,3],:
）=[]%删除指定行
D1=
04040404
04040404
【例3.2-8】函数flipud,fliplr,rot90对数组的操作体现着“矩阵变换”。
A=reshape（1:
9,3,3）
A=
147
258
369
B=flipud（A）%上下对称交换
B=
369
258
147
C=fliplr（A）%左右对称交换
C=
741
852
963
D=rot90（A,2）%逆时针旋转90度，2次
D=
963
852
741
.3数组、矩阵及向量化编程
MATLAB面向数组/矩阵编程和运算：
Ø用“数组或矩阵运算”模式去处理那些“借助循环而反复执行的标量运算”
●显著提高程序执行速度
●书写简洁、便于阅读
3.3.1数组概念和运算规则
（3-2）
数组运算符和矩阵运算符
加
减
乘
左除
右除
求幂
数组运算符
.+
.-
.*
.\
./
.^
矩阵运算符
+
-
*
\
/
^
服从数组运算规则的MATLAB初等函数及关系逻辑算符
举例
服从数组运算规则的
初等函数
三角、反三角
sin,cos,tan,cot,sec,csc,asin,acos,atan,acos,asec,acsc
双曲、反双曲
sinh,cosh,…,asinh,acosh,…
指数、对数
exp,sqrt,pow2,log,log10,log2
圆整、求余
ceil,floor,fix,round,mod,rem
模、角、虚实部
abs,angle,real,imag,conj
符号函数
sign
关系运算符
==,~=,>,<,>=,<=
逻辑运算符
&,|,~
3.3.2矩阵概念及矩阵运算规则
10一矩阵概念
（3-1）
10二矩阵运算规则
【例】已知矩阵
，
，采用两种不同的编程求这两个矩阵的乘积
。
（1）以“标量”为基本处理单元的传统编程法
clear
rng（'default'）%随机发生器恢复默认设置
A=rand（2,4）;
B=rand（4,3）;
C1=zeros（size（A,1）,size（B,2））;
forii=1:
size（A,1）%<6>
forjj=1:
size（B,2）
fork=1:
size（A,2）
C1（ii,jj）=C1（ii,jj）+A（ii,k）*B（k,jj）;
end
end
end%<12>
C1
C1=
1.27261.38701.2281
2.29481.46601.8204
（2）MATLAB求矩阵乘积
C2=A*B%<13>
C2=
1.27261.38701.2281
2.29481.46601.8204
3.3.3数组、矩阵算术运算规则的功能比较
相同：
加，减，与标量加、减、乘
不同：
转置，乘，除，求幂
【例】本例演示：
矩阵和数组乘法的不同。
（1）
Am=magic（3）%3阶魔方阵
Aa=reshape（1:
12,3,4）
B=2*ones（3,4）
Am=
816
357
492
Aa=
14710
25811
36912
B=
2222
2222
2222
（2）
Cm1=Am*B
Ca1=Aa.*B
Cm1=
30303030
30303030
30303030
Ca1=
281420
4101622
6121824
3.3.4向量化编程
尽可能用“数组或矩阵运算”指令
【例】欧姆定律：
，其中
分别是电阻（欧姆）、电压（伏特）、电流（安培）。
验证实验：
据电阻两端施加的电压，测量电阻中流过的电流，然后据测得的电压、电流计算平均电阻值。
（测得的电压电流具体数据见下列程序）。
（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）;
fork=1:
L
r（k）=vr（k）/ir（k）;
end
%---------------------------
sr=0;
fork=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）%MATLAB现成的求平均函数
r=
Columns1through7
31.785730.000030.900029.448330.678029.883730.0669
Columns8through10
30.817131.493530.1739
rm=
30.5247
【例】用间距为0.1的水平线和垂直线均匀分割
的矩形域，在所有水平线和垂直线交点上计算函数
的值，并图示。
（1）非向量化编程
clear
x=-5:
0.1:
5;
y=（-2.5:
0.1:
2.5）';
N=length（x）;%101
M=length（y）;%51
forii=1:
M%51行
forjj=1:
N%101列
X0（ii,jj）=x（jj）;%所有格点的x坐标
Y0（ii,jj）=y（ii）;%所有格点的y坐标
Z0（ii,jj）=sin（abs（x（jj）*y（ii）））;%所有格点的函数值
end
end
（2）向量化编程
[X,Y]=meshgrid（x,y）;%指定矩形域内所有格点的（x,y）坐标
Z=sin（abs（X.*Y））;%数组运算计算矩形域所有格点坐标（x,y）对应的函数值
%注意：
函数f（·）对数组的逐个元素起作用。
（3）比较二维双精度数是否相等
norm（Z-Z0）%范数接近eps,认为相等。
ans=
0
（4）绘图
surf（X,Y,Z）
xlabel（'x'）
ylabel（'y'）
shadinginterp
view（[190,70]）
图3.3-1指定域上的二元函数图形
surf（X0,Y0,Z0）
xlabel（'x'）
ylabel（'y'）
shadinginterp
view（[190,70]）
3.4“非数”和“空”数组
✧MATLAB中特有的两个概念和“预定义变量”
3.4.1非数NaN（或记为nan）
由
，
，
，∞-∞等运算产生。
NaN的性质：
●NaN参与运算所得的结果也是NaN，即具有传递性；
●NaN没有“大小”概念，不能比较两个NaN的大小。
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）%将随机发生器置0
R=rand（2,5）;R（1,5）=NaN;R（2,3）=NaN
R=
0.95010.60680.89130.4565NaN
0.23110.4860NaN0.01850.4447
LR=isnan（R）%对数组元素是否非数进行判断
LR=
00001
00100
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个元素
补充：
find函数寻找真值1，并给出相应的坐标位置。
Ind2sub（size（a）,i）函数给出在a数组中单下标为i的数值的全下标。
相反，sub2ind（size（a）,x,y）给出在a数组中坐标为（x,y）的数值的单下标。
3.4.2“空”数组
Ø定义：
数组的某维长度为0或若干维长度均为0。
●MATLAB为操作和表述需要专门设计的一种数组。
【例3.4-3】关于“空”数组的算例。
（1）创建“空”数组的几种方法
a=[]%二维“空”（行、列长度均为0）数组
b=ones（2,0）
c=zeros（2,0）%“空”数组不同于全零数组
d=eye（2,0）
f=rand（2,3,0,4）
a=
[]
b=
Emptymatrix:
2-by-0
c=
Emptymatrix:
2-by-0
d=
Emptymatrix:
2-by-0
f=
Emptyarray:
2-by-3-by-0-by-4
（2）“空”数组的属性
class（a）%数据类型
isnumeric（a）%是数值数组吗？
isempty（a）%唯一可正确判断数组是否“空”的指令
ans=
double
ans=
1
ans=
1
whicha%a是什么
ndims（a）%维数
size（a）%数组大小注意：
“空”数组并非“虚无”，它确实存在。
aisavariable.
ans=
2
ans=
00
补充：
size（a）是给出数组a有几行几列，而ndims（a）则给出数组a有几维。
（3）“空”数组用于子数组的删除和大数组的大小收缩
A=reshape（-4:
5,2,5）
A=
-4-2024
-3-1135
A（:
[2,4]）=[]
A=
-404
-315
3.5关系操作和逻辑操作
MATLAB约定：
●关系（或逻辑）表达式中，作为输入的任何非0数都被看作是“逻辑真”，而只有0才被认为是“逻辑假”。
●计算结果，即输出，是由0和1组成的“逻辑数组（LogicalArray）”。
在此数组中的1表示“真”，0表示“假”。
●逻辑数组是一种特殊的数值数组。
与“数值类”有关的操作和函数
升级会员 