第二讲 Matlab的数值计算 —— Matlab 具有出色的数值计算能力,占据世界上数值计算软件的主导地位。 —— Matlab数值计算是使用Matlab的基础,是Matlab强大计算功能的体现
一、Matlab的表达式和变量 1、Matlab的表达式 Matlab语句最常见的两种形式: 表达式 变量=表达式
对于第二种形式,系统将右边的表达式的结果赋值给左边的变量,然后在屏幕上显示出来 对于第一种形式,计算结果系统自动赋值给名为ans的变量,ans永远只保存最近一次的表达式的运算结果 对于第二种形式,系统将右边的表达式的结果赋值给左边的变量,然后在屏幕上显示出来 注意:表达式末尾的;
数字表达 Matlab的数值采用十进制表示,可以使用科学计数法,用e表示位数。 Matlab常用运算符 + - * /(右除) \(左除)^(幂)
2、Matlab的变量 Matlab对变量名的大小写敏感,大小写代表不同的变量. 变量名可以是字母或数字,但首字必须是字母. Matlab对变量名的大小写敏感,大小写代表不同的变量. 显示结果的缺省格式是5位有效数字,可以用命令format改变输出格式,如:format short e; format long e。 Matlab可以响应键盘输入,用“input”命令: x=input(‘please input x: ’)
二、Matlab的基本计算功能 常用基本数学函数 函数名称 函数功能 abs(x) 取绝对值 sign(x) 符号函数 angle(z) 复数的相角 rem(x,y) x/y取余 sqrt(x) 开平方 gcd(x,y) 最大公因数 real(z) 复数的实部 lcm(x,y) 最小公倍数 imag(z) 复数的虚部 exp(x) 自然指数 conj(z) 共轭复数 pow2(x) 2的指数 round(x) 四舍五入取整 log(x) 自然对数 flx(x) 无论正负,舍去小数取整 log2(x) 以2为底的对数 rat(x) 将实数化为分数 log10(x) 以10为底的对数
有:sin(x),cos (x),tan (x), asin (x),acos (x),atan (x)等 Matlab常用的三角函数 有:sin(x),cos (x),tan (x), asin (x),acos (x),atan (x)等
三、矩阵的创建 1、创建矩阵的方法 直接输入法 规则: 矩阵元素必须用[ ]括住 矩阵元素必须用逗号或空格分隔 矩阵元素必须用[ ]括住 矩阵元素必须用逗号或空格分隔 在[ ]内矩阵的行与行之间必须 用分号分隔
矩阵元素 a=[1 2 3;4 5 6] x=[2 pi/2;sqrt(3) 3+5i] 矩阵元素可以是任何Matlab表达式 ,可以是实数 ,也可以是复数,复数可用特殊函数i,j 输入 a=[1 2 3;4 5 6] x=[2 pi/2;sqrt(3) 3+5i]
符号的作用 逗号和分号的作用 逗号和分号可作为指令间的分隔符,Matlab允许多条语句在同一行出现。 分号如果出现在指令后,屏幕上将不显示结果。
注意:只要是赋过值的变量,不管是否在屏幕上显示过,都存储在工作空间中,以后可随时显示或调用。变量名尽可能不要重复,否则会覆盖 。 当一个指令或矩阵太长时,可用•••续行
冒号的作用 用于生成等间隔的向量,默认间隔为1。 用于选出矩阵指定行、列及元素。 循环语句
2.用Matlab函数创建矩阵 空阵 [ ] — Matlab允许输入空阵,当一项操作无结果时,返回空阵。 rand —— 随机矩阵 eye —— 单位矩阵 zeros ——全部元素都为0的矩阵 ones ——全部元素都为1的矩阵 diag ——产生对角阵
例:0~1分布的随机矩阵,用rand函数可以产生任意行列的的0~1分布的随机矩阵 >> ra=rand(2,3) ra = 0.9218 0.1763 0.9355 0.7382 0.4057 0.9169
注意:Matlab严格区分大小写字母,因此a与A是两个不同的变量。 还有伴随矩阵、稀疏矩阵、魔方矩阵、对角矩阵、范德蒙等矩阵的创建,就不一一介绍了。 注意:Matlab严格区分大小写字母,因此a与A是两个不同的变量。 Matlab函数名必须小写。
3. 矩阵的修改 直接修改 可用键找到所要修改的矩阵,用键移动到要修改的矩阵元素上即可修改。 指令修改 直接修改 可用键找到所要修改的矩阵,用键移动到要修改的矩阵元素上即可修改。 指令修改 可以用A(,)= 来修改。
3)由矩阵编辑器修改 由Matlab提供工具栏按钮来查看工作区变量 单击变量,可以打开或 删除变量 打开后得到如下图, 点击矩阵元素,修改即可
例如 a=[1 2 0;3 0 5;7 8 9] a =1 2 0 3 0 5 7 8 9 a(3,3)=0 7 8 0
四、数据的保存与获取 把Matlab工作空间中一些有用的数据长久保存下来的方法是生成mat数据文件。 save —— 将工作空间中所有的变量存到Matlab.mat文件中。
save data——将工作空间中所有的变量存到data.mat文件中。 save data a b ——将工作空间中a和b变量存到data.mat文件中。 下次运行Matlab时即可用load指令调用已生成的mat文件。
mat文件是标准的二进制文件,还可以ASCII码形式保存。 load —— load data —— load data a b —— mat文件是标准的二进制文件,还可以ASCII码形式保存。 即可恢复保存过的所有变量
五、矩阵运算 1.矩阵加、减(+,-)运算 规则: 相加、减的两矩阵必须有相同的行和列两矩阵对应元素相加减。 允许参与运算的两矩阵之一是标量。标量与矩阵的所有元素分别进行加减操作。
2. 矩阵乘()运算 规则: A矩阵的列数必须等于B矩阵的行数 标量可与任何矩阵相乘。 a=[1 2 3;4 5 6;7 8 0];b=[1;2;3];c=a*b c =14 32 23
d=[-1;0;2];f=pi*d f = -3.1416 6.2832
3. 矩阵除( \,/)运算 矩阵除的运算在线性代数中没有,有矩阵逆的运算,在Matlab中有两种矩阵除运算: \左除 /右除 A\B为方程AX=B的解 B/A为方程XA=B的解 A\B=inv(A)*B 而B/A=B*inv(A) B/A=(A’\B’)’
其中inv(B)表示矩阵B的逆。
4. 矩阵乘方—— a^n,a^p,p^a a ^ p —— a 自乘p次幂 对于p的其它值,计算将涉及特征值 方阵 >1的整数 对于p的其它值,计算将涉及特征值 和特征向量,如果p是矩阵,a是标量 a^p使用特征值和特征向量自乘到p次 幂;如a,p都是矩阵,a^p则无意义。
※当一个方阵有复数特征值或负实特征值时,非整数幂是复数阵。 a=[1,2,3;4,5,6;7,8,9];a^2 ans =30 36 42 66 81 96 102 126 150 ※当一个方阵有复数特征值或负实特征值时,非整数幂是复数阵。
a^0.5 ans = 0.4498 + 0.7623i 0.5526 + 0.2068i 0.6555 -0.3487i 1.0185 + 0.0842i 1.2515 + 0.0228i 1.4844 - 0.0385i 1.5873 - 0.5940i 1.9503 - 0.1611i 2.3134 + 0.2717i
5. 矩阵的其它运算 ’ —— 矩阵转置(共轭)
6.矩阵标识和子矩阵 矩阵的子矩阵可以通过标量、向量、冒号的标识来引用和赋值。 子阵序号向量标识方式A(v,w) 矩阵的子矩阵可以通过标量、向量、冒号的标识来引用和赋值。 子阵序号向量标识方式A(v,w) v,w可以是任何排列的向量,也可以是“:”,它表示全部行(在V的位置)或全部列(在w的位置)。 “0~1”向量标识方式A(L1,:)、 A(:,L2)、 A(L1,L2) 向量L1,L2的长度分别为矩阵A的行数和列数,向量L1,L2中的元素取1(表示提取相应的行或列)或0(不提取)
例:提取矩阵的子阵 B=magic(5) B1=B(1:2,[1,3,5]) B2=B([3,1],:) B([1,3],[2,4])=zeros(2)
6.矩阵的一些特殊操作 矩阵的变维 a=[1:12],b=reshape(a,3,4) %b为3行4列 c=zeros(3,4);c(:)=a(:) 矩阵的变向 rot90:旋转(逆时针); fliplr:左右翻; flipud:上下翻 矩阵的抽取 diag:抽取主对角线;tril: 抽取主下三角; triu:抽取主上三角
例:a=[1 2 3;4 5 6;7 8 9] , b=rot90(a),c=flipud(a),d=fliplr(a) a = 1 2 3 4 5 6 7 8 9 b = 3 6 9 2 5 8 1 4 7 c = 7 8 9 4 5 6 1 2 3 d = 3 2 1 6 5 4 9 8 7
la=tril(a,-1),ua=triu(a,1),na=triu(a) 4 0 0 7 8 0 ua = 0 2 3 0 0 6 0 0 0 na = 1 2 3 0 5 6 0 0 9
7、特殊的应用矩阵 在控制系统分析和设计中,应用一些特殊的矩阵 约当阵--将状态空间模型转换为约当标准形 形式:jordan(a) 得到矩阵a的约当标准形 [V,J]=jordan(a) J为约当标准形,V为相似变换阵,满足V\a*V=J. 例:a=[1,2;3,4];J=jordan(a),[V,JJ]=jordan(a) J = 5.3723 0 0 -0.3723 V = 0.2389 0.7611 0.5222 -0.5222 JJ = 5.3723 0
伴随矩阵--将状态空间模型转换为可控标准形 compan(p) 获得多项式P的伴随矩阵 p=[1 2 3],a=compan(p) p = 1 2 3 a = -2 -3 1 0
格雷姆(Gram)矩阵-计算可控性和可观性格雷姆阵 Gc=gram(G,’c’) 计算稳定系统G的可控性格雷姆矩阵Gc Go=gram(G,’o’) 计算稳定系统G的观测性性格雷姆矩阵Go
7.常用的矩阵函数 矩阵的行列式、矩阵的秩、特征值等在现代控制理论中有广泛的应用,Matlab提供了相应的函数求其值 det(A) 方阵A的行列式 eig(A) 方阵A的特征值和特征向量 rank(A) 矩阵A的秩 trace(A) 矩阵A的迹 expm(A) 矩阵的指数 sqrtm(A) 求矩阵的平方根 funm(A,’fun’) 求一般的方阵函数
A=[1 2 3;4 5 6;7 8 0] deta=det(A) ranka=rank(A) tracea=trace(A) eiga=eig(A) >>deta =27 ranka =3 tracea = 6 eiga = 12.1229 -0.3884 -5.7345
注意:funm函数的问题 f=funm(A,’fun’)对于方阵A,计算用函数‘fun’定义的矩阵函数,例如:funm(A,’sin’) 就是计算矩阵A的sin函数。 funm使用的算法有潜在的不稳定性,如果A是接近于一个多重特殊值和病态特征向量的矩阵,它会产生不精确的值,Matlab系统会因此出现警告的信息,但是有时候出现正确结果时,也会出现警告信息。
例:矩阵函数和数组函数运算的区别 B=[1/6,1/2;1/3,5/6]*pi >> sinb=sin(B) sinb = 0.5000 1.0000 0.8660 0.5000 >> fsinmb=funm(B,'sin') fsinmb = -0.0849 0.0000 0.0000 -0.0849
8. 矩阵的数组运算 与通常意义上的由符号表示的线性代数 矩阵运算不同 在Matlab中提供了如下的数组运算符: 数组运算指元素对元素的算术运算, 与通常意义上的由符号表示的线性代数 矩阵运算不同 在Matlab中提供了如下的数组运算符: .+ .- .* .^ .\(./) .‘
注意:运算符中的小黑点绝对不能遗漏,否则将不按数组运算规则进行计算。不管执行什么数组计算,所计算结果数组总是与参与运算的数组同维。数组运算中所有的二元运算必须是同维的数组或者其中有一个是标量。 1. 数组加减(.+,.-) a.+b a.- b
2. 数组乘除(,./,.\) ab —— a,b两数组必须有相同的行 和列两数组相应元素相乘。 ans = 2 8 18 4 15 30 49 72 90
a=[1 2 3;4 5 6;7 8 9]; b=[2 4 6;1 3 5;7 9 10]; a*b ans = 25 37 46 25 37 46 55 85 109 85 133 172
a./b=b.\a — 都是a的元素被b的对应元 素除 a.\b=b./a — 都是b的元素被a的对应元 例: a=[1 2 3];b=[4 5 6]; c1=a.\b; c2=b./a c1 = 4.0000 2.5000 2.0000 c2 = 4.0000 2.5000 2.0000 —— 给出a,b对应元素间的商.
3. 数组乘方(.^) — 元素对元素的幂 例: a=[1 2 3];b=[4 5 6]; z=a.^2 z = 3. 数组乘方(.^) — 元素对元素的幂 例: a=[1 2 3];b=[4 5 6]; z=a.^2 z = 1.00 4.00 9.00 z=a.^b 1.00 32.00 729.00
4. 数组转置(.‘) — 元素对元素的转置 >>a=[1+2*i 3*i;4+5*i 6*i],aa=a.',aaa=a' a = 1.0000 + 2.0000i 0 + 3.0000i 4.0000 + 5.0000i 0 + 6.0000i aa = 1.0000 + 2.0000i 4.0000 + 5.0000i 0 + 3.0000i 0 + 6.0000i aaa = 1.0000 - 2.0000i 4.0000 - 5.0000i 0 - 3.0000i 0 - 6.0000i
六、 多项式运算 Matlab语言把多项式表达成一个行向量, 该向量中的元素是按多项式降幂排列的。 f(x)=anxn+an-1xn-1+……+a0 可用行向量 p=[an an-1 …… a1 a0]表示 1. poly —— 产生特征多项式系数向量 特征多项式一定是n+1维的 特征多项式第一个元素一定是1
例:a=[1 2 3;4 5 6;7 8 0]; p=poly(a) p =1.00 -6.00 -72.00 -27.00 p是多项式p(x)=x3-6x2-72x-27的Matlab描述方法,我们可用: p1=poly2str(p,'x') — 函数文件,显示 数学多项式的形式 p1 =x^3 - 6 x^2 - 72 x - 27
2.roots —— 求多项式的根 a=[1 2 3;4 5 6;7 8 0];p=poly(a) p = 1.00 -6.00 -72.00 -27.00 r=roots(p) r = 12.12 -5.73 ——显然 r是矩阵a的特征值 -0.39
当然我们可用poly令其返回多项式形式 p2=poly(r) p2 = 1.00 -6.00 -72.00 -27.00 Matlab规定多项式系数向量用行向量表示,一组根用列向量表示
3.conv,convs多项式乘运算 例:a(x)=x2+2x+3; b(x)=4x2+5x+6; c = (x2+2x+3)(4x2+5x+6) a=[1 2 3];b=[4 5 6]; c=conv(a,b)=conv([1 2 3],[4 5 6]) c = 4.00 13.00 28.00 27.00 18.00 p=poly2str(c,'x') p = 4 x^4 + 13 x^3 + 28 x^2 + 27 x + 18
4.deconv多项式除运算 a=[1 2 3]; c = [4.00 13.00 28.00 27.00 18.00] d=deconv(c,a) d =4.00 5.00 6.00 [d,r]=deconv(c,a) 余数 c除a后的整数
5.多项式微分 Matlab提供了polyder函数多项式的微分。 命令格式: polyder(p): 求p的微分 polyder(a,b): 求多项式a,b乘积的微分 [p,q]=polyder(a,b): 求多项式a,b商的微分,即有理多项式求导 例:a=[1 2 3 4 5]; poly2str(a,'x') ans = x^4 + 2 x^3 + 3 x^2 + 4 x + 5 b=polyder(a) b = 4 6 6 4 poly2str(b,'x') ans =4 x^3 + 6 x^2 + 6 x + 4
例:有理多项式求导 num=10*[1 4 5 6 7] den=poly([-2;-1;-0.5;]) [p,q]=polyder(num,den) %num分子多项式, den分母多项式 p = 10 70 195 200 -125 -390 -185 q = 1.0000 7.0000 19.2500 26.5000 19.2500 7.0000 1.0000 该结果证实:
Matlab提供了函数residue执行部分分式展开。 6.多项式展开与组合 Matlab提供了函数residue执行部分分式展开。 命令格式: [r,p,k]=residue(b,a)或者相应的逆运算[b,a]=residue[r,p,k] 例: num=10*[1 4 5 6 7] den=poly([-2;-1;-0.5;]) %根的形式 [res,poles,k]=residue(num,den) [n,d]=residue(res,poles,k)
7.多项式估值 按数组运算规则计算多项式的值 PA=polyval(p,S) p为多项式,S为矩阵,矩阵S中每一个值代入多项式求值,结果是与S同维的矩阵。 按矩阵规则计算多项式的值 PAM=polyvalm(p,S) p为多项式,S为矩阵,矩阵S作为变量代入多项式求出值,结果是与S同维的矩阵。
例:多项式计算 s=rand(3),p=poly(s), pa=polyval(p,s),pm=polyvalm(p,s) s = 0.9501 0.4860 0.4565 0.2311 0.8913 0.0185 0.6068 0.7621 0.8214 p = 1.0000 -2.6628 1.9560 -0.4289 pa =-0.1166 0.0075 0.0042 -0.1067 -0.0929 -0.3936 0.0009 -0.0422 -0.0647 pm = 1.0e-015 * -0.3331 0.1665 -0.0555 -0.0260 -0.0555 -0.0451 -0.1665 0.2776 -0.3331
s=rand(3),a=1;p=poly(s),pa=polyval(p,a), pm=polyvalm(p,a) s = 0.6038 0.0153 0.9318 0.2722 0.7468 0.4660 0.1988 0.4451 0.4186 p = 1.0000 -1.7692 0.6195 -0.0378 pa = -0.1875 pm =