《MATLAB的数值计算1.ppt》由会员分享，可在线阅读，更多相关《MATLAB的数值计算1.ppt（69页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、一、命令行的基本操作创建矩阵的方法创建矩阵的方法直接输入法规则： 矩阵元素必须用 括住 矩阵元素必须用逗号或空格分隔 在 内矩阵的行与行之间必须 用分号分隔 矩阵元素可以是任何matlab表达式 ，可以是实数 ，也可以是复数，复数可用特殊函数I，j 输入 a=1 2 3;4 5 6 x=2 pi/2;sqrt(3) 3+5i 矩阵元素符号的作用逗号和分号的作用 逗号和分号可作为指令间的分隔符，matlab允许多条语句在同一行出现。 分号如果出现在指令后，屏幕上将不显示结果。注意：只要是赋过值的变量，不管是否在屏幕上显示过，都存储在工作空间中，以后可随时显示或调用。变量名尽可能不要重复，否则会覆
2、盖 。 当一个指令或矩阵太长时，可用续行冒号的作用 用于生成等间隔的向量，默认间隔为1。 用于选出矩阵指定行、列及元素。 循环语句2.用matlab函数创建矩阵空阵 matlab允许输入空阵，当一项操作无结果时，返回空阵。rand 随机矩阵eye 单位矩阵zeros 全部元素都为0的矩阵ones 全部元素都为1的矩阵 还有伴随矩阵、稀疏矩阵、魔方矩阵、对角矩阵、范德蒙等矩阵的创建，就不一一介绍了。注意：matlab严格区分大小写字母，因此a与A是两个不同的变量。 matlab函数名必须小写。3. 矩阵的修改 直接修改 可用键找到所要修改的矩阵，用键移动到要修改的矩阵元素上即可修改。 指令修改
3、可以用A(,)= 来修改。例如例如a=1 2 0;3 0 5;7 8 9a =1 2 0 3 0 5 7 8 9a(3,3)=0a =1 2 0 3 0 5 7 8 0把matlab工作空间中一些有用的数据长久保存下来的方法是生成mat数据文件。 save 将工作空间中所有的变量存到matlab.mat文件中。二、数据的保存与获取默认文件名save data将工作空间中所有的变量存到data.mat文件中。save data a b 将工作空间中a和b变量存到data.mat文件中。 下次运行matlab时即可用load指令调用已生成的mat文件。load load data load dat
4、a a b mat文件是标准的二进制文件，还可以ASCII码形式保存。即可恢复保存过的所有变量矩阵加、减（,）运算规则： 相加、减的两矩阵必须有相同的行和列两矩阵对应元素相加减。 允许参与运算的两矩阵之一是标量。标量与矩阵的所有元素分别进行加减操作。三、矩阵运算2. 矩阵乘（）运算规则：A矩阵的列数必须等于B矩阵的行数标量可与任何矩阵相乘。a=1 2 3;4 5 6;7 8 0;b=1;2;3;c=a*bc =14 32 23 d=-1;0;2;f=pi*df = -3.1416 0 6.2832 矩阵除的运算在线性代数中没有，有矩阵逆的运算，在matlab中有两种矩阵除运算 a p a 自乘
5、p次幂 方阵方阵1的整数的整数3. 矩阵乘方 an,ap,pa对于p的其它值,计算将涉及特征值和特征向量，如果p是矩阵，a是标量ap使用特征值和特征向量自乘到p次幂；如a,p都是矩阵，ap则无意义。 a=1,2,3;4,5,6;7,8,9;a2 ans =30 36 42 66 81 96 102 126 150当一个方阵有复数特征值或负实特征值时，非整数幂是复数阵。 a0.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.
6、5873 - 0.5940i 1.9503 - 0.1611i 2.3134 + 0.2717iinv 矩阵求逆det 行列式的值eig 矩阵的特征值diag 对角矩阵 矩阵转置sqrt 矩阵开方4. 矩阵的其它运算 5.矩阵的一些特殊操作矩阵的变维 a=1:12;b=reshape(a,3,4) c=zeros(3,4);c(:)=a(:)矩阵的变向 rot90:旋转; fliplr:上翻; flipud:下翻矩阵的抽取 diag:抽取主对角线;tril: 抽取主下三角; triu:抽取主上三角矩阵的扩展关系运算 关系符号意义=小于小于或等于大于大于或等于等于不等于 数组运算指元素对元素的算
7、术运算，与通常意义上的由符号表示的线性代数矩阵运算不同 数组加减(.+,.-) a.+b a.- b5. 矩阵的数组运算 对应元素相加减（与矩阵加对应元素相加减（与矩阵加减等效）减等效）2. 数组乘除(，./，.）ab a，b两数组必须有相同的行 和列两数组相应元素相乘。a=1 2 3;4 5 6;7 8 9;b=2 4 6;1 3 5;7 9 10;a.*bans = 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*bans = 25 37 46 55 85 109 85 133 172 a./b=b.aa
8、.b=b./aa./b=b.a 都是a的元素被b的对应元 素除a.b=b./a 都是a的元素被b的对应元 素除例: a=1 2 3;b=4 5 6; c1=a.b; c2=b./ac1 = 4.0000 2.5000 2.0000c2 = 4.0000 2.5000 2.0000 给出a,b对应元素间的商.3. 数组乘方(.) 元素对元素的幂例:a=1 2 3;b=4 5 6;z=a.2z = 1.00 4.00 9.00z=a.bz = 1.00 32.00 729.00matlab语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。 f(x)=anxn+an-1xn-1+lo
9、a0 可用行向量 p=an an-1 a1 +a0表示poly 产生特征多项式系数向量特征多项式一定是n+1维的1.特征多项式第一个元素一定是1四、 多项式运算 在MATLAB中，多项式使用降幂系数的行向量表示，如：多项式11625012234xxxxp=poly(r)p = 1 -12 -0 25 116（1）多项式的建立与表示方法r=roots(p)r = 11.7473 2.7028 -1.2251 + 1.4672i -1.2251 - 1.4672i表示为：p=1 -12 0 25 116，使用函数roots可以求出多项式等于0的根，根用列向量表示。若已知多项式等于0的根，函数pol
10、y可以求出相应多项式。例: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 =x3 - 6 x2 - 72 x - 27（2）多项式的运算相乘conva=1 2 3 ; b=1 2 c=conv(a,b)=1 4 7 6conv指令可以嵌套使用，如conv(conv(a,b),c)相除deconvq,r=deconv(c,b)q=1 2 3 商多项式r=0 0 0 余多项式求多项式的
11、微分多项式polyderpolyder(a)=2 2求多项式函数值polyval(p,n)：将值n代入多项式求解。polyval(a,2)=11conv,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.00p=poly2str(c,x)p = 4 x4 + 13 x3 + 28 x2 + 27 x + 18deconv多项式除运算a=1 2 3; c = 4.00 1
12、3.00 28.00 27.00 18.00d=deconv(c,a)d =4.00 5.00 6.00d,r=deconv(c,a)余数余数c除除a后的整数后的整数多项式微分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 = x4 + 2 x3 + 3 x2 + 4 x + 5b=polyder(a)b = 4 6 6 4poly2str(b,x)ans =4 x3 + 6
13、x2 + 6 x + 4（3）*多项式的拟合多项式拟合又称为曲线拟合，其目的就是在众多的样本点中进行拟合，找出满足样本点分布的多项式。这在分析实验数据，将实验数据做解析描述时非常有用。命令格式：p=polyfit(x,y,n)，其中x和y为样本点向量，n为所求多项式的阶数，p为求出的多项式。例exp2_15.m%curve fitting of sin waveclcclearx=0:0.1:2*pi; %生成样本点xy=sin(x)+0.5*rand(size(x); %生成样本点y，通过随机矩阵p=polyfit(x,y,3) %拟合出多项式（3阶）y1=polyval(p,x); %求多
14、项式的值plot(x,y,+,x,y1,-r) %绘制多项式曲线，以验证结果4）*多项式插值多项式插值是指根据给定的有限个样本点，产生另外的估计点以达到数据更为平滑的效果。该技巧在信号处理与图像处理上应用广泛。所用指令有一维的interp1、二维的interp2、三维的interp3。这些指令分别有不同的方法（method），设计者可以根据需要选择适当的方法，以满足系统属性的要求。Help polyfun可以得到更详细的内容。y=interp1(xs,ys,x,method)在有限样本点向量xs与ys中，插值产生向量x和y，所用方法定义在method中，有4种选择：nearest：执行速度最快
15、，输出结果为直角转折linear：默认值，在样本点上斜率变化很大spline：最花时间，但输出结果也最平滑cubic：最占内存，输出结果与spline差不多例exp2_16.m%curve interpolationys=0 0.9 0.6 1 0 0.1 -0.3 -0.7 -0.9 -0.2; %已有的样本点ysxs=0:length(ys)-1; %已有的样本点xsx=0:0.1:length(ys)-1;%新的样本点xy1=interp1(xs,ys,x,nearest); %插值产生新的样本点y1y2=interp1(xs,ys,x,linear); %插值产生新的样本点y2y3=i
16、nterp1(xs,ys,x,spline); %插值产生新的样本点y3y4=interp1(xs,ys,x,cubic); %插值产生新的样本点y4plot(xs,ys,+k,x,y1,:r,x,y2,-m,x,y3,-c,x,y4,-b) %分别绘制不同方法产生的曲线legend(sampled point,nearest,linear,spline,cubic)五、代数方程组求解matlab中有两种除运算左除和右除。对于方程ax+b，a 为anm矩阵，有三种情况： 当n=m时，此方程成为“恰定”方程 当nm时，此方程成为“超定”方程 当nm时，此方程成为“欠定”方程 matlab定义的除
17、运算可以很方便地解上述三种方程1.恰定方程组的解方程ax+b(a为非奇异) x=a-1 b 矩阵逆两种解:x=inv(a)b 采用求逆运算解方程 x=ab 采用左除运算解方程 方程ax=ba=1 2;2 3;b=8;13;x=inv(a)*b x=ab x = x = 2.00 2.00 3.00 3.00322121xx138 = a x = b例: x1+2x2=8 2x1+3x2=132.超定方程组的解方程 ax=b ,mn时此时不存在唯一解。方程解 (a a)x=a b x=(a a)-1 a b 求逆法 x=ab matlab用最小二乘法找一 个准确地基本解。 例: x1+2x2=1
18、 2x1+3x2=2 3x1+4x2=3a=1 2;2 3;3 4;b=1;2;3; 解1 x=ab 解2 x=inv(aa) a b x = x = 1.00 1.00 0 0.00 21xx321 =433221 a x = b3.欠定方程组的解 当方程数少于未知量个数时,即不定情况,有无穷多个解存在。matlab可求出两个解：用除法求的解x是具有最多零元素的解是具有最小长度或范数的解，这个解是基于伪逆pinv求得的。 x1+2x2+3x3=1 2x1+3x2+4x3=2a=1 2 3;2 3 4;b=1;2; x=ab x=pinv(a)b x = x = 1.00 0.83 0 0.3
19、3 0 -0.17432321321xxx21=a x = b4、矩阵分解（1）奇异值分解U,S,V=svd(A)例：a = 9 8 6 8可以验证：u*u=Iv*v=Iu*s*v=a求矩阵A的奇异值及分解矩阵，满足U*S*V=A，其中U、V矩阵为正交矩阵（U*U=I），S矩阵为对角矩阵，它的对角元素即A矩阵的奇异值。u,s,v=svd(a)u = 0.7705 -0.6375 0.6375 0.7705s = 15.5765 0 0 1.5408v = 0.6907 -0.7231 0.7231 0.6907（2）特征值分解V,D=eig(A)例： a = 9 8 6 8v,d=eig(a)
20、v = 0.7787 -0.7320 0.6274 0.6813d = 15.4462 0 0 1.5538求矩阵A的特征向量V及特征值D，满足A*V=V*D。其中D的对角线元素为特征值，V的列为对应的特征向量。如果D=eig(A)则只返回特征值。可以验证：A*V=V*D（3）正交分解Q,R=qr(A)例： a = 9 8 6 8q,r=qr(a)q = -0.8321 -0.5547 -0.5547 0.8321r = -10.8167 -11.0940 0 2.2188将矩阵A做正交化分解，使得Q*R=A，其中Q为正交矩阵（其范数为1，指令norm(Q)=1)，R为对角化的上三角矩阵。no
21、rm(q) ans = 1q*rans = 9.0000 8.0000 6.0000 8.0000（4）三角分解L,U=lu(A)将A做对角线分解，使得A=L*U,其中L为下三角矩阵，U为上三角矩阵。注意：L实际上是一个“心理上”的下三角矩阵，它事实上是一个置换矩阵P的逆矩阵与一个真正下三角矩阵L1（其对角线元素为1）的乘积。L1,U1,P=lu(A)例：a=1 2 3;4 5 6;7 8 9 比较： l1,u1,p=lu(a) l,u=lu(a)l1 = 1.00 0 0 0.14 1.00 0 0.57 0.50 1.00u1 = 7.00 8.00 9.00 0 0.86 1.71 0
22、0 0.00p = 0 0 1 1 0 0 0 1 0l = 0.14 1.00 0 0.57 0.50 1.00 1.00 0 0u = 7.00 8.00 9.00 0 0.86 1.71 0 0 0.00可以验证：u1=u，inv(p)*l1=la=l*up*a=l1*u1 数值积分法数值积分法 一般情况下，在控制系统仿真中最常用、最基本的求解常微分方程数值解的方法主要是数值积分法。 对于形如 的系统，已知系统状态变量y的初值y0,现要计算y随时间变化的过程y(t)，对微分方程的积分可以写作：右图所示曲线下的面积就是y(t),由于难以得到积分的数值表达式，所以采用近似的方法，常用有三种形
23、式：欧拉法梯形法龙格一库塔法 ttdtytftyy0),()(0),(tuyfy 六、微分方程求解),()(11kkkkkythfyyyt欧拉公式欧拉公式,采用矩形面积近似积分结果，即当t= tk+1时hk=tk+1-tk,若步距不变，则hk=h.为了提高精度，人们提出了“梯形法”，其中最简单的是亚当姆斯预报亚当姆斯预报校正公式校正公式 ，先用欧拉法估计近似值，然后用梯形法进行校正：),(),(21),(1)0(111)0(kkkkkkkkkkytfytfhyyythfyy龙格-库塔法的基本原理在连续系统仿真中，主要的数值计算工作是求解一阶微分方程：),( ytfdtdy若已知y的初值y0,再
24、按离散化原理，对上式我们可以写成：ttkkdtytftytykk1),()()(1再对上式的右端函数f(t,y)(为任意非线性函数)在tk附近展开成泰勒级数，依照展开的阶次不同我们就构成了不同的龙格-库塔公式。二阶龙格二阶龙格库塔公式库塔公式,记在tk时刻的状态变量为yk,并定义两个附加向量型变量 :),( ),( )(2121211hkyhtfkytfkkkhyykkkkkk四阶龙格库塔公式 ：),()2,2()2,2(),()22(6342312143211hkyhtfkkhyhtfkkhyhtfkytfkkkkkhyykkkkkkkkkk不论几阶RK法，它们的计算公式都是由两部分组成，即
25、上一步的结果yk和步长h乖以tk至tk+1时间间隔间各外推点的导数ki的加权平均和有了上面的数学算法，就可以用MATLA编写出该算法的函数： functiontout,yout = rk4(odefile,tspan,y0) t0= tspan(1);th= tspan(2); if length(tspan)= 3,h= tspan(3) else, h= tspan(2)-tspan(1);th= tspan(2); end tout= t0:h:th;yout= ; for t= tout k1= eval(odefile(t,y0); k2= eval(odefile(t+h/2,y0
26、+0.5*k1); k3= eval(odefile(t+h/2,y0+0.5*k2); k4= eval(odefilet+h,y0+k3); yout= yout;y0;end其中，tspan可以有两种构成方法：一是可以是一个等间距的时间向量；二是tspan=t0,th,h,t0和th为计算的初始及终止值，h为计算步长，odefile是一个字符串变量，表示微分方程的文件名，y0是初值列向量，函数调用完成后，时间向量与各个时刻状态变量构成的矩阵分别由tout和yout返回. MATLAB提供了两个RK法函数ode23()和ode45()，一个采用二阶三级公式，一个采用四阶五级RK法，并采用了
27、自适应变步长的求解方法，即当解的变化较慢时采用较大的计算步长，以加快计算速度；当方程的解变化得较快时，积分步长自动变小，以得到较高的计算精度。例：x+(x2-1)x+x=0为方便令x1=x，x2=x分别对x1,x2求一阶导数，整理后写成一阶微分方程组形式 x1=x2 x2=x2(1-x12)-x1建立m文件1.解微分方程建立m文件function xdot=wf(t,x)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,
28、 x0)plot(t,x), figure(2),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)-2.5-2-1.5-1-0.500.511.522.5-3-2-10123051015202530-3-2-10123七、函数优化寻优函数：fmin 单变量函数fmins 多变量函数constr
29、有约束条件无约束条件无约束条件例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)八、数据分析max 各列最大值各列最大值 mean 各列平均值各列平均值sum 各列求和各列求和std 各列标准差各列标准差var 各列方差各列方差sort 各列递增排序各列递增排序小 结 本节介绍了matlab语言的数值运算功能
30、，通过学习应该掌握：如何创建矩阵、修改矩阵符号的用法矩阵及数组运算多项式运算线性方程组与微分运算进入夏天，少不了一个热字当头，电扇空调陆续登场，每逢此时，总会进入夏天，少不了一个热字当头，电扇空调陆续登场，每逢此时，总会想起那一把蒲扇。蒲扇，是记忆中的农村，夏季经常用的一件物品。记想起那一把蒲扇。蒲扇，是记忆中的农村，夏季经常用的一件物品。记忆中的故乡，每逢进入夏天，集市上最常见的便是蒲扇、凉席，不论男女老忆中的故乡，每逢进入夏天，集市上最常见的便是蒲扇、凉席，不论男女老少，个个手持一把，忽闪忽闪个不停，嘴里叨叨着少，个个手持一把，忽闪忽闪个不停，嘴里叨叨着“怎么这么热怎么这么热”，于是三，于
31、是三五成群，聚在大树下，或站着，或随即坐在石头上，手持那把扇子，边唠嗑五成群，聚在大树下，或站着，或随即坐在石头上，手持那把扇子，边唠嗑边乘凉。孩子们却在周围跑跑跳跳，热得满头大汗，不时听到边乘凉。孩子们却在周围跑跑跳跳，热得满头大汗，不时听到“强子，别跑强子，别跑了，快来我给你扇扇了，快来我给你扇扇”。孩子们才不听这一套，跑个没完，直到累气喘吁吁，。孩子们才不听这一套，跑个没完，直到累气喘吁吁，这才一跑一踮地围过了，这时母亲总是，好似生气的样子，边扇边训，这才一跑一踮地围过了，这时母亲总是，好似生气的样子，边扇边训，“你你看热的，跑什么？看热的，跑什么？”此时这把蒲扇，是那么凉快，那么的温馨
32、幸福，有母亲此时这把蒲扇，是那么凉快，那么的温馨幸福，有母亲的味道！蒲扇是中国传统工艺品，在我国已有三千年多年的历史。取材的味道！蒲扇是中国传统工艺品，在我国已有三千年多年的历史。取材于棕榈树，制作简单，方便携带，且蒲扇的表面光滑，因而，古人常会在上于棕榈树，制作简单，方便携带，且蒲扇的表面光滑，因而，古人常会在上面作画。古有棕扇、葵扇、蒲扇、蕉扇诸名，实即今日的蒲扇，江浙称之为面作画。古有棕扇、葵扇、蒲扇、蕉扇诸名，实即今日的蒲扇，江浙称之为芭蕉扇。六七十年代，人们最常用的就是这种，似圆非圆，轻巧又便宜的蒲芭蕉扇。六七十年代，人们最常用的就是这种，似圆非圆，轻巧又便宜的蒲扇。蒲扇流传至今，我的记忆中，它跨越了半个世纪，也走过了我们的扇。蒲扇流传至今，我的记忆中，它跨越了半个世纪，也走过了我们的半个人生的轨迹，携带着特有的念想，一年年，一天天，流向长长的时间隧半个人生的轨迹，携带着特有的念想，一年年，一天天，流向长长的时间隧道，袅道，袅