MATLAB的数值运算.ppt

1、第二章 MATLAB的数值运算,第一节 基本语法结构,1.变量,命名,字母+任意字母(数字、下划线),规则,字母的大小写、标点符号,变量操作,命令窗口,命令、变量值,存储,调用,一、变量与赋值,2赋值语句,用运算符将运算量连接 结果是矩阵,变量=表达式,赋值,回车,语句执行,显示执行结果,; 结果不显示,二、预定义变量,工作空间中，系统本身定义。 Pi、i，j 特定的含义，避免重新赋值。,z = -0.3488 + 0.3286i,例 计算表达式的值，并显示计算结果。,x=1+2i;y=3-sqrt(17);z=(cos(abs(x+y)-sin(78*pi/180)/(x+abs(y),屏幕

2、显示，赋值给指定变量 没有指定变量时，赋值给特殊变量ans 显示格式- format Format-结果的显示 MATLAB总是以双字长浮点数(双精度)来执行所有的运算。,三、数值显示格式,结果为整数，显示没有小数 结果不是整数，输出形式： format (short)：短格式（5位定点数）99.1253 format long：长格式（15位定点数 99.12345678900000) format short e：短格式e方式 9.9123e+001 format long e：长格式e方式 9.912345678900000e+001 format bank：2位十进制 99.12 fo

3、rmat hex：十六进制格式,三、数值显示格式,四、 常用数学函数,使用说明 三角函数以弧度为单位计算 abs函数-求实数的绝对值、复数的模、字符串的ASCII码值,自变量,数学函数,运算法则,运算结果,矩阵变量,函数逐项作用于矩阵的元素,与自变量同维数的矩阵,第二节矩阵运算,一、创建矩阵的方法,1、直接输入法,y=2,4,5;3 6 8 y= 2 4 5 3 6 8,a=1; b=2; c=3; x=5 b c; a*b a+c c/b x= 5.000 2.000 3.000 2.000 4.000 1.500, 逗号或空格 分号,任何matlab表达式 实数 ，复数 复数可用特殊函数i

4、，j 输入,矩阵元素,a=1 2 3;4 5 6 x=2 pi/2;sqrt(3) 3+5i,2.函数创建矩阵,5 0 0 0 7 0 0 0 2,1 0 0 0 1 0,0 0 0 0 0 0,1 1 1 1 1 1,eye(2,3) zeros(2,3) ones(2,3) V=5 7 2; A=diag(V),1 0 0 1 0 0 0 0 1 1 1 1,eye(2) zeros(2) ones(2),单位矩阵 零矩阵 一矩阵 对角矩阵,空阵 :当操作无结果时，返回空阵。 对角矩阵：对角元素向量 V=a1,a2,an A=diag(V) A为方阵，V=diag(A)提取A的对角元素构成

5、向量V。 随机矩阵：rand(m,n) mn的均匀分布,3利用冒号表达式建立向量,a=1:2:10,初始值 ：步长： 终止值,e1: e2: e3,a= 1 3 5 7 9,第一个元素，最后元素，元素总数,linspace(a,b,n),=a:(b-a)/(n-1):b,a=linspace(1,10,10),a= 1 2 3 4 5 6 7 8 9 10,a(3,3)=0,4. 矩阵的修改,a=1 2 0;3 0 5;7 8 9 a =1 2 0 3 0 5 7 8 9,直接修改- ,指令修改-A(,)= ,a =1 2 0 3 0 5 7 8 0,a = 1 4 2 5 3 6,a = 1

6、 4 2 5 3 6,二、矩阵的运算,1.转置, .,a=1 2 3;4 5 6,a=1 2 3;4 5 6.,实矩阵,转置结果相同,复数矩阵,b = 1.0000 + 2.0000i 2.0000 - 7.0000i,b=1+2i 2-7i,b=1+2i 2-7i.,b = 1.0000 - 2.0000i 2.0000 + 7.0000i,对复数进行共轭处理 将其排列形式进行转置,标量矩阵，标量与矩阵所有元素分别加减操作,2.矩阵加、减（,）运算,a=1 2 0;3 0 5;7 8 9;,b=1 2 0;3 0 5;7 8 9;,c=1;,相同的行和列对应元素相加减,a+b;a-b;,a+

7、c;a-c;,a=1 2 3;4 5 6;7 8 0 ; b=1;2;3;c=a*b,A矩阵的列数必须等于B矩阵的行数 标量可与任何矩阵相乘,d=-1;0;2; f=pi*d,c =14 32 23,f = -3.1416 0 6.2832,3.矩阵乘（）运算,a p a 自乘p次幂,方阵,1的整数,如果p是矩阵，a是标量： ap使用特征值和特征向量自乘到p次幂； 如a,p都是矩阵，则ap无意义。,4.矩阵乘方 an,ap,pa,a=1,2,3;4,5,6;7,8,9; a2 ans =30 36 42 66 81 96 102 126 150,5.矩阵除法,inv 矩阵求逆 det 行列式的

8、值 eig 矩阵的特征值 diag 对角矩阵 sqrt 矩阵开方,6.矩阵的其它运算,在有关算术运算符前面加点-.*、./、.、.,矩阵的对应元素进行相关运算 要求两矩阵的维参数相同,7. 矩阵的点(数组)运算,a=1 2 3;4 5 6;7 8 9; b=2 4 6;1 3 5;7 9 10;,a.*b a*b,例: 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=b.a a.b=b./a a./b=b.a a.b=b./a, 给出a,b对应元素间的

9、商., 都是a的元素被b的对应元素除, 都是a的元素被b的对应元素除,a=1 2 3;b=4 5 6; z1=a.2；z2=a.b,z1 = 1.00 4.00 9.00,z 2= 1.00 32.00 729.00,数组乘方(.) 元素对元素的幂,求逆：inv(A)； 求行列式：det(A),-2.3333 0.3333 1.0000 2.6667 0.3333 -2.0000 -0.6667 -0.3333 1.0000,9.逆矩阵与行列式计算,矩阵必须为方阵,a=1 2 3; 4 5 6; 2 3 5; b=inv(a) det(a),-3,a=1 2;3 4；b= 3 5; 5 9 c

10、=a+b；d=a-b a*b=13 23; 29 51 a/b=-0.50 0.50;3.50 1.50 ab=-1 -1;2 3 a3=37 54; 81 118 a.*b=3 10;15 36 a./b=0.33 0.40;0.60 0.44 a.b=3.00 2.50;1.67 2.25 a.3= 1 8; 27 64,例题,三、矩阵的操作,1、矩阵下标,提取子块,矩阵扩展,重排子块,A(m,n),A(:,n),A(m,:),A(m1:m2,n1:n2),A(:), ,在矩阵中一个不存在的地址位置上赋值,m = 2 n = 3,ans = 3 max(size(a) ans = 3,an

11、s = 2,2、矩阵的大小,a=1 2 3;3 4 5; m,n=size(a),length(a),rank(a),矩阵的行列数,行数或列数的最大值,矩阵的秩,3、矩阵分解,奇异值分解,U*S*V=A,正交矩阵,对角矩阵,U,S,V=svd(A),a =9 8;6 8; u,s,v=svd(a),u = 0.7705 -0.6375 0.6375 0.7705,s = 15.5765 0 0 1.5408 v = 0.6907 -0.7231 0.7231 0.6907,u*u=I v*v=I u*s*v=a,特征值分解,a =9 8;6 8; v,d=eig(a),若D=eig(A)，只返

12、回特征值,A*V=V*D,A的特征向量,A的特征值,列为特征值,对角线元素为特征向量,V,D=eig(A),v = 0.7787 -0.7320 0.6274 0.6813 d = 15.4462 0 0 1.5538,正交三角形分解,a=9 8;6 8; q,r=qr(a),Q,R=qr(A),q = -0.8321 -0.5547 -0.5547 0.8321 r = -10.8167 -11.0940 0 2.2188,正交矩阵,上三角矩阵,上三角分解,L,U=lu(A),置换矩阵,上三角矩阵,将方阵A做上三角分解，使得A=L*U,a=1 2 3;4 5 6;7 8 9; l,u=lu(

13、a),l = 0.14 1.00 0 0.57 0.50 1.00 1.00 0 0,u = 7.00 8.00 9.00 0 0.86 1.71 0 0 0.00,多项式,元素按多项式降幂排列,f(x)=anxn+an-1xn-1+a0,一、多项式的建立与表示方法,第三节 多项式运算,一个行向量,p=an an-1 a1 a0,已知多项式等于0的根，求出相应多项式,r = 11.7473 2.7028 -1.2251 + 1.4672i -1.2251 - 1.4672i,p = 1 -12 -0 25 116,p=1 -12 0 25 116,roots,poly,多项式等于0的根，列向量

14、,r=roots(p),p=poly(r),nomial, 函数文件，显示数学多项式,特征多项式-n+1维 特征多项式-第一个元素是1,a=1 2 3;4 5 6;7 8 0; p=poly(a),p1 =x3 - 6 x2 - 72 x - 27,1、poly 产生特征多项式系数向量,p =1.00 -6.00 -72.00 -27.00,p(x)=x3-6x2-72x-27,p1=poly2str(p,x),p =1.00 -6.00 -72.00 -27.00,r = 12.12 -5.73 -0.39,返回多项式形式,p2=poly(r) p2 =1.00 -6.00 -72.00 -

15、27.00,2、roots 求多项式的根,a=1 2 3;4 5 6;7 8 0; p=poly(a),r=roots(p),r是矩阵a的特征值,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 x4 + 13 x3 + 28 x2 + 27 x + 18,1、多项式乘运算,conv指令嵌套使用,二、多项式运算,conv(P1,P2),多项式系数向量,conv(conv(a,b),c),a(x)=x2+2x+3;b(x)=4x2+5x+6;,c = (x2+2x+3)(4x2

16、+5x+6),a=1 2 3;b=4 5 6;,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(p1,p2),p1除以p2的商式,deconv是conv的逆函数，P1=conv(P2,d)+r。,2、deconv多项式除运算,余数,多项式系数向量,polyder(p) polyder(a,b) p,q=polyder(a,b),a=1 2 3 4 5; poly2str(a,x) ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyder(a) b

17、 = 4 6 6 4 poly2str(b,x) ans =4 x3 + 6 x2 + 6 x + 4,3、多项式微分,求p的微分,求多项式a,b乘积的微分,求多项式a,b商的微分,第四节代数方程组求解,matlab定义的除运算(左除、右除)可以很方便地解上述三种方程,方程ax=b,“恰定”方程,“超定”方程,“欠定”方程,n=m,nm,nm,anm矩阵,一、恰定方程组的解,方程ax=b(a为非奇异),x=a-1 b,矩阵逆,x=inv(a)b 求逆运算解方程,x=ab 左除运算解方程,=,a x = b,x1+2x2=8 2x1+3x2=13,1 2 2 3,x1 x2,8 13,a=1 2

18、;2 3;b=8;13; x=inv(a)*b x=ab,x = 2.00 3.00,二、超定方程组的解,方程ax=b,anm矩阵,mn时不存在唯一解,(a a)x=a b,x=(a a)-1 a b 求逆法,x=ab 最小二乘法准确基本解,x1+2x2=1 2x1+3x2=2 3x1+4x2=3,=,a x = b,1 2 3 3 4,x1 x2,1 2 3,a=1 2;2 3;3 4;b=1;2;3;, x=ab x=inv(aa) ab,x = 1.00 0.00,三、欠定方程组的解,方程ax=b,anm矩阵,mn,存在无穷多解,方程数少于未知量个数,用除法求的解x,基于伪逆pinv求得

19、,具有最多零元素,具有最小长度或范数,x1+2x2+3x3=1 2x1+3x2+4x3=2,1 2 3 2 3 4,x1 X2 x3,=,1 2,a x = b,a=1 2 3;2 3 4;b=1;2;, x=ab x=pinv(a)b, x = 1.00 0.00 0.00, x =0.83 0.33 -0.17,第五节 统计分析,一、最大值和最小值max、min,1向量的最大值、最小值,y=max(X),向量X的最大值，若X中有复数元素，按模取最大值,x=-43,72,9,16,23,47; y=max(x)； y,l=max(x),y,I=max(X),最大值的序号,y=min(X),y

20、,I=min(X),2求矩阵的最大值和最小值,max(A),Y,U=max(A),行向量，列最大值,行向量，列最大值,行向量，列最大值的行号,max(A,dim),dim=1,同max(A),dim=2,列向量，行大值,3两个向量或矩阵对应元素的比较,U=max(A,B),U=max(A,n),同型向量或矩阵,与AB同型向量或矩阵， A,B对应元素的较大者,标量,与A同型向量或矩阵， A对应元素和n中较大者,U=min(A,B),U=min(A,n),二、求和与求积,X-向量，A-矩阵,sum(X),-X各元素和,prod(X),-积,sum(A),-行向量，各列元素的和,prod(A),-行

21、向量，各列元素的积,sum(A,dim),dim=1,同sum(A),dim=2,列向量，各行元素和,prod(A,dim),dim=1,同prod(A),dim=2,列向量，各行元素积,三、 平均值和中值,mean(X),median(X),X-向量，A-矩阵,-X的算术平均值,-X的中值,mean(A)-行向量，各列的算术平均值 median(A)-行向量，各列的中值,四、累加和与累乘积,cumsum(X)-向量X累加和向量 cumprod(X)-向量X累乘积向量 cumsum(A)-矩阵，第i列是A的第i列累加和向量 cumprod(A)-矩阵，第i列是A的第i列累乘积向量,五、标准方差

22、与相关系数,缺省flag=0，dim=1,1求标准方差,std,Y=std(A,flag,dim),2相关系数,例：生成满足正态分布的100005随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。,X=randn(10000,5); M=mean(X) D=std(X) R=corrcoef(X),corrcoef(X)-从矩阵X形成的一个相关系数矩阵。大小同矩阵X。它把矩阵X的每列作为一个变量，然后求其相关系数。,Y,I=sort(A,dim)-对各列或各行重新排序,Y中的元素在A中位置,六、 排序,sort(X)-按升序排列的新向量,排序后的矩阵,定义对某些集合给

23、定的数据点之间函数的估值。 求出所需中间点的函数 一维、二维、 三次样条等许多插值,第六节 数据插值,插值函数 table1、table2、interp1、interp2 、spline,X1的取值范围不能超出X的给定范围， 否则，会给出“NaN”错误。,一、一维数据插值,Y1=interp1(X,Y,X1,method),两个等长已知向量，采样点、样本值,Y1=interp1(X,Y,X1,method),插值方法 linear、nearest、cubic、spline,插值结果 与X1等长,欲插值的点，向量或标量,设时间变量h为一行向量，温度变量t为一个两列矩阵，第一列存放室内温度，第二列

24、储存室外温度。 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),Y1=spline(X,Y,X1)- 3次样条插值函数,例：某观测站测得某日6:00时至18:00时之间每隔2小时的室内外温度()，用3次样条插值分别求得该日室内外6:30至17:30时之间每隔2小时各点的近似温度()。,二、二维数据插值,X1、Y1的取值范围不能超出X的给定范围， 否则，会给出“NaN”错误。,Z1=interp2(X,Y,Z,X1,Y1,method),向量或矩阵，采

25、样点,插值方法 linear、nearest、cubic、spline,插值结果,欲插值的点，向量或标量,与参数采样点对应的函数值,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),某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度()。试用线性插值求出在1分钟内每隔20秒、钢轨每隔1米处的温度TI。,fmin 单变量函数

26、 fmins 多变量函数 constr 有约束条件,无约束条件,第七节函数优化,例1：f(x)=x2+3x+2在-5 5区间的最小值 f=fmin(x2+3*x+2,-5,5) 例2：f(x)=100(x2-x12)2+(a-x1)2在x1=a, x2=a2处有最小值 function f=xun(x,a) f=100*(x(2)-x(1).2).2+(a-x(1).2; x=fmins(xun,0,0, , ,sqrt(2),第八节微分方程求解,Euler（欧拉法）、Runge Kutta(龙格-库塔法）。 Euler法-一步法，用于一阶微分方程,所以 yn+1 = yn + hf (xn,

27、yn) n=0,1,2 y(x0)=y0,当给定仿真步长时：,龙格-库塔法-实际上取两点斜率的平均斜率来计算的，其精度高于欧拉算法 。 龙格-库塔法：ode23 ode45,Runge Kutta法,例：x+(x2-1)x+x=0,建立m文件 解微分方程,function xdot=wf(t,x)%建立m文件 xdot=zeros(2,1) xdot(1)=x(2) xdot(2)=x(2)*(1-x(1)2)-x(1) 给定区间、初始值;求解微分方程 t0=0; tf=20; x0=0 0.25; t,x=ode23(wf, t0, tf, x0) plot(t,x), figure(2),

28、plot(x(:,1),x(:,2),T,Y = ODE23(ODEFUN,TSPAN,Y0) 建立m文件 function dxdt=wf(t,x) dxdt=x(2);x(2)*(1-x(1)2)-x(1); 求解微分方程 t,x=ode23(wf,0 30,0 0.25); plot(t,x); figure(2) plot(x(:,1),x(:,2),第九节拟合,x0=0:0.1:1； y0=-.447 1.978 3.11 5.25 5.02 4.66 4.01 4.58 3.45 5.35 9.22; p=polyfit(x0,y0,3) p = 56.6915 -87.1174

29、40.0070 -0.9043 xx=0:0.01:1;yy=polyval(p,xx); plot(xx,yy,-b,x0,y0,or),多项式拟合,(1) 建立矩阵A。A=4,-65,-54,0,6;56,0,67,-45,0(2) 找出大于4的元素的位置。find(A4),1. 建立矩阵A，然后找出大于4的元素的位置。,习题,2. 先建立55矩阵A，然后将A的第一行元素乘以1，第二行乘以2，第五行乘以5。,A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3; 11,18,25,2,19;D=diag(1:5);D*A %用D左乘A，对

30、A的每行乘以一个指定常数,3x5-7x4+5x2+2x-18=0p=3,-7,0,5,2,-18;A=compan(p); %A的伴随矩阵x1=eig(A) %求A的特征值x2=roots(p) %直接求多项式p的零点,3 用求特征值的方法解方程。,4 分别建立33、32和与矩阵A同样大小的零矩阵。 (1) 建立一个33零矩阵。zeros(3) (2) 建立一个32零矩阵。zeros(3,2) (3) 设A为23矩阵，则可以用zeros(size(A)建立一个与矩阵A同样大小零矩阵。A=1 2 3;4 5 6; %产生一个23阶矩阵Azeros(size(A) %产生一个与矩阵A同样大小的零矩阵,5 求多项式x4+8x3-10与多项式2x2-x+3的乘积。 6 求多项式