1、第二讲 MATLAB的数值计算, matlab 具有出色的数值计算能力，占据世界上数值计算软件的主导地位,数值运算的功能,创建矩阵 矩阵运算 多项式运算 线性方程组 数值统计 线性插值 函数优化 微分方程的数值解,一、命令行的基本操作,创建矩阵的方法 直接输入法 规则： 矩阵元素必须用 括住 矩阵元素必须用逗号或空格分隔 在 内矩阵的行与行之间必须 用分号分隔,矩阵元素可以是任何matlab表达式 ，可以是实数 ，也可以是复数，复数可用特殊函数i，j 输入 a=1 2 3;4 5 6x=2 pi/2;sqrt(3) 3+5i,矩阵元素,符号的作用,逗号和分号的作用逗号和分号可作为指令间的分隔符
2、，matlab允许多条语句在同一行出现。分号如果出现在指令后，屏幕上将不显示结果。,注意：只要是赋过值的变量，不管是否在屏幕上显示过，都存储在工作空间中，以后可随时显示或调用。变量名尽可能不要重复，否则会覆盖 。当一个指令或矩阵太长时，可用续行,冒号的作用用于生成等间隔的向量，默认间隔为1。用于选出矩阵指定行、列及元素。循环语句,2.用matlab函数创建矩阵,空阵 matlab允许输入空阵，当一项操作无结果时，返回空阵。 rand 随机矩阵 eye 单位矩阵 zeros 全部元素都为0的矩阵 ones 全部元素都为1的矩阵,还有伴随矩阵、稀疏矩阵、魔方矩阵、对角矩阵、范德蒙等矩阵的创建，就不
3、一一介绍了。 注意：matlab严格区分大小写字母，因此a与A是两个不同的变量。matlab函数名必须小写。,3. 矩阵的修改, 直接修改可用键找到所要修改的矩阵，用键移动到要修改的矩阵元素上即可修改。 指令修改可以用A(,)= 来修改。,例如 a=1 2 0;3 0 5;7 8 9 a =1 2 03 0 57 8 9 a(3,3)=0 a =1 2 03 0 57 8 0,把matlab工作空间中一些有用的数据长久保存下来的方法是生成mat数据文件。save 将工作空间中所有的变量存到matlab.mat文件中。,二、数据的保存与获取,默认文件名,save data将工作空间中所有的变量存
4、到data.mat文件中。 save data a b 将工作空间中a和b变量存到data.mat文件中。下次运行matlab时即可用load指令调用已生成的mat文件。,load load data load data a b mat文件是标准的二进制文件，还可以ASCII码形式保存。,即可恢复保存过的所有变量,矩阵加、减（,）运算 规则： 相加、减的两矩阵必须有相同的行和列两矩阵对应元素相加减。 允许参与运算的两矩阵之一是标量。标量与矩阵的所有元素分别进行加减操作。,三、矩阵运算,如： A=1,2;3,4; b=3Ab =4 56 7,2. 矩阵乘（）运算 规则： A矩阵的列数必须等于B矩
5、阵的行数 标量可与任何矩阵相乘。 a=1 2 3;4 5 6;7 8 0;b=1;2;3;c=a*b c =143223,d=-1;0;2;f=pi*df = -3.141606.2832矩阵除的运算在线性代数中没有，有矩阵逆的运算，在matlab中有两种矩阵除运算, a=1,2;3,4;b=1,0;0,1;a/b %a右除b ans =1 2 a*inv(b)3 4 a=1,2;3,4;b=1,0;0,1;ab a左除b ans =-2.0000 1.0000 % inv (a)*b1.5000 -0.5000 ,a p a 自乘p次幂,方阵,1的整数,3. 矩阵乘方 an,ap,pa,对于
6、p的其它值,计算将涉及特征值 和特征向量，如果p是矩阵，a是标量 ap使用特征值和特征向量自乘到p次 幂；如a,p都是矩阵，ap则无意义。,a=1,2,3;4,5,6;7,8,9;a2ans =30 36 4266 81 96102 126 150,当一个方阵有复数特征值或负实特征值时，非整数幂是复数阵。,a0.5ans =0.4498 + 0.7623i 0.5526 + 0.2068i 0.6555 -0.3487i1.0185 + 0.0842i 1.2515 + 0.0228i 1.4844 - 0.0385i1.5873 - 0.5940i 1.9503 - 0.1611i 2.31
7、34 + 0.2717i,inv 矩阵求逆 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:抽取主上三角； 矩阵的扩展a=1,2;3,4a(6,6)=1,关系运算, a=1:9 a =1 2 3 4 5 6 7 8 9 b=10-a b =9 8 7 6 5 4 3
8、 2 1 r0=(a r1=(a=b) r1 =0 0 0 0 1 0 0 0 0,数组运算指元素对元素的算术运算， 与通常意义上的由符号表示的线性代数 矩阵运算不同数组加减(.+,.-)a.+ba.- b,5. 矩阵的数组运算,对应元素相加减（与矩阵加减等效）,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.*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 a
9、ns =25 37 46 55 85 109 85 133 172,a./b=b.a a.b=b./a a./b=b.a 都是a的元素被b的对应元 素除 a除b a.b=b./a 都是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 =1.00 4.00 9.00 z=a.b =14 25 36 z =1.00 3
10、2.00 729.0014 25 36,matlab语言把多项式表达成一个行向量， 该向量中的元素是按多项式降幂排列的。f(x)=anxn+an-1xn-1+a0可用行向量 p=an an-1 a1 a0表示,四、 多项式运算,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.00p是多项式p(x)=x3-6x2-72x-27的matlab描述方法，我们可用： p1=poly2str(p,x) 函数文件，显示 数学多项式的形式 p1 =x
11、3 - 6 x2 - 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，多项式乘运算,例:a(x)=x2+2x+3; b(x)=4x2+5x+6; c = (x2+2x+3)(4x2+5x+6)
12、a=1 2 3;b=4 5 6; c=conv(a,b)或c=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,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,5.多项式微分,matlab提供了polyder函数多项式的微分。 命令格式： polyder(p): 求p的微分 polyder(a,b): 求多项式a,b乘积的微分 p,q=po
13、lyder(a,b): 求多项式a,b商的微分，p/q,例：a=1 2 3 4 5; poly2str(a,x) ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyder(a) b = 4 6 6 4 poly2str(b,x) ans =4 x3 + 6 x2 + 6 x + 4,五、代数方程组求解,matlab中有两种除运算左除和右除。 对于方程ax=b，a 为anm矩阵，有三种情况： 当n=m时，此方程成为“恰定”方程 当nm时，此方程成为“超定”方程 当nm时，此方程成为“欠定”方程matlab定义的除运算可以很方便地解上 述三种方程,1.恰定方程组的解,方
14、程ax=b(a为非奇异)x=a-1 b 矩阵逆 两种解: x=inv(a)b 采用求逆运算解方程 x=ab 采用左除运算解方程,方程ax=b a=1 2;2 3;b=8;13; x=inv(a)*b x=abx = x = 2.00 2.003.00 3.00,=,a x = b,例: x1+2x2=8 2x1+3x2=13,2.超定方程组的解,方程 ax=b ,mn时此时不存在唯一解。 方程解 (a a)x=a b x=(a a)-1 a b 求逆法 x=ab matlab用最小二乘法找一个准确地基本解。,例: x1+2x2=1 2x1+3x2=23x1+4x2=3a=1 2;2 3;3 4
15、;b=1;2;3;解1 x=ab 解2 x=inv(aa) a b x = x =1.00 1.000 0.00,=,a x = b,3.欠定方程组的解,当方程数少于未知量个数时,即不定 情况,有无穷多个解存在。 matlab可求出两个解： 用除法求的解x是具有最多零元素的解。 是具有最小长度或范数的解，这个解是基于伪逆pinv求得的。,x1+2x2+3x3=12x1+3x2+4x3=2a=1 2 3;2 3 4;b=1;2;x=ab x=pinv(a)b x = x =1.00 0.830 0.330 -0.17,a x = b,六、微分方程求解,微分方程求解的仿真算法有多种，常用的有Eul
16、er（欧拉法）、Runge Kutta(龙格-库塔法。 Euler法称一步法，用于一阶微分方程,当给定仿真步长时：所以 yn+1 = yn + hf (xn,yn) n=0,1,2y(x0)=y0,Runge Kutta法龙格-库塔法：实际上取两点斜率 的平均 斜率来计算的，其精度高 于欧拉算法 。 龙格-库塔法：ode23 ode45,k1=hf(xn,yn) k2=hf(xn+h,yn+k),例：x+(x2-1)x+x=0 为方便令x1=x，x2=x分别对x1,x2求一 阶导数，整理后写成一阶微分方程组 形式x1=x2x2=x2(1-x12)-x1 建立m文件 解微分方程,M文件,当待解决
17、问题的指令较多且所用指令较复杂时，或当一组指令通过改变少量参数就可以被反复使用去解决不同问题时。建立M脚本文件建立运行,M文件及程序调试与变量属性,M文件：由MATLAB命令语句构成的程序文件称作M文件，它将m作为文件的扩展名。可以直接阅读并用任何文本编辑器来建立这种文件。 M文件可分为两种：主程序和子程序。主程序 (也称主程序文件)： 一般是由用户为解决特定的问题而编制的程序。子程序 (称为函数文件)：它必须由其它M文件来调用。函数文件往往具有一定的通用性，并且可以进行递归调用。,M文件的结构,(1)有关程序的功能、使用方法等内容的注释。主程序前面的若干行通常是程序的注释， 每行以“”开始。
18、 (2)程序的主体。 (3)整个程序应按MATLAB标识符的要求起文件名，并加上后缀m。,对主程序文件的运行方式通常有两种：,在MATLAB的命令窗口中运行：在MATLAB的命令窗中键入此程序的文件名后，系统就开始执行文件中的程序； 在编辑窗口中运行：通过编辑窗口打开所要运行的文件，然后再运行。主程序文件中的语句可以对MATLAB工作空间中的所有数据进行运算操作。,函数文件结构,函数文件是用来定义子程序的。它与主程序文件的主要区别有三点： (1)由function起头，后跟的函数名必须与文件名相同； (2)有输入输出变元(变量)，可进行变量传递； (3)除非用global声明，程序中的变量均为
19、局部变量，不保存在工作空间中。 通常，函数文件由五部分构成： 函数定义行； H1行； 函数帮助文本； 函数体； 注释。,函数定义命令function( ),格式： function output1,output2,=functionname(input1,input2,.)功能：定义函数名为functionname的函数(或子程序)，该函数的输出变量为output1,output2,，输出变量为 (input1,input2,.)。,例如：function y=mean(x,dim) function为函数定义的关键字，mean为函数名，y为输出变量，x和dim为输入变量。 注意：当函数具有多
20、个输出变量时，则以方括号括起；当函数具有多个输入变量时，则直接用圆括号括起。例如：functionx，y，zsphere(theta，phi，rho)。当函数不含输出变量时，则直接略去输出部分或采用空方括号表示。,程序的调试,1.MATLAB的主程序的调试：MATLAB的查错能力很强，加上工作空间中变量的保存和显示功能不需要专门的调试命令，因此，主程序的调试很方便地进行。,2.函数程序的调试：,因为在函数程序中出错而停机时，其变量不作保存。虽然它也会指出出错的语句，但因为子程序中的变量(局部变量)在程序执行完毕后会自动消失，其它现场数据又无记录，会给调试带来很大困难。可采用下列措施： (1)把
21、某些分号改为逗号，使中间结果能显示在屏幕上，作为查错的依据。 (2)在子程序中适当部位加keyboard命令，到了此处，系统会暂停而等待用户键入命令。这时子程序中的变量还存于工作空间中，可以对它进行检查。 (3)将函数文件的第一行前加“%”号，使它成为程序文件来作初步调试。第一行中的输入变元，可改用input或赋值语句来输入，调试好后再改回为函数文件。 (4)在程序中使用以下人机交互命令，也可大大提高程序的调试程度。,建立m文件 function xdot=wf(t,x) xdot=zeros(2,1) xdot(1)=x(2) xdot(2)=x(2)*(1-x(1)2)-x(1) 输在文件
22、中 给定区间、初始值;求解微分方程 t0=0; tf=20; x0=0 0.25; t,x=ode23(wf, t0, tf, 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),七、函数优化 寻优函数： fminbnd 单变量函
23、数 fminsearch 多变量函数 constr 有约束条件,无约束条件,例1：f(x)=x2+3x+2在-5 5区间的最小值 f=fminbnd(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=fminsearch(xun,0,0, , ,sqrt(2),八、数据分析与插值函数,max 各列最大值 mean 各列平均值 sum 各列求和 std 各列标准差 var 各列方差 sort 各列递增排序,九、拟合与
24、插值,1. 多项式拟合 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 40.0070 -0.9043 xx=0:0.01:1;yy=polyval(p,xx); plot(xx,yy,-b,x0,y0,or),2.插值插值的定义是对某些集合给定的数据点之间函数的估值方法。 当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具。 Matlab提供了一维、二维、 三次样条等许多插值选择,table1 table
25、2 intep1 interp2 spline 利用已知点确定未知点 粗糙 精确 集合大的 简化的,插值函数,小 结本节介绍了matlab语言的数值运算 功能，通过学习应该掌握： 如何创建矩阵、修改矩阵 符号的用法 矩阵及数组运算 多项式运算 线性方程组与微分运算,基本数学函数,三角函数1正弦sin函数与反正弦函数asin函数 格式：Asin(X) 功能：对矩阵或数组X的每个元素求正弦值，所有角度用弧度表示。 格式：Xasin(A) 功能：对于A的元素求反正弦值，结果为弧度。,2双曲正弦sinh函数与反双曲正弦函数asinh函数 格式：Asinh(X) 功能：对矩阵或数组X的每个元素求双曲正弦
26、值。 格式：Xasinh(A) 功能：对于A的元素求反双曲正弦值。,3余弦cos函数与反余弦函数acos函数 格式：Acos(X) 功能：对矩阵或数组X的每个元素求余弦值，所有角度用弧度表示。 格式：Xacos(A) 功能：对于A的元素求反余弦值，结果为弧度。,4双曲余弦cosh函数与反双曲余弦函数acosh函数 格式：Acosh(X) 功能：对矩阵或数组X的元素求双曲余弦值。 格式：Xacosh(A) 功能：对于A的元素求反双曲余弦值。,5正切tan函数与反正切函数atan函数 格式：Atan(X) 功能：对矩阵或数组X的每个元素求正切值，所有角度用弧度表示。 格式：Xatan(A) 功能：
27、对于A的每个元素求反正切值，结果为弧度。,6四象限反正切atan2函数 格式：P=atan2(Y，X) 功能：四象限反正切函数。返回大小与Y、X相同的矩阵或数组，由Y和X对应元素的 实部求反正切值得到，其中虚部忽略。,7双曲正切tanh函数与反双曲正切atanh函数 格式：Atanh(X) 功能：对于X的每个元素求双曲正切值。 格式：Xatanh(A) 功能：对于A的每个元素求反双曲正切值。,8正割sec函数与反正割asec函数 格式：Asec(X) 功能：对X的每个元素求正割。 格式：Xasec(A) 功能：返回A的每个元素的反正割值。,9双曲正割sec函数与反双曲正割asec函数 格式：A
28、sech(X) 功能：对X的每个元素求双曲正割。 格式：Xasech(A) 功能：返回A的每个元素的反双曲正割值。,10余割csc函数与反余割acsc函数 格式：Acsc(X) 功能：对X的每个元素求正割。 格式：Xachc(A) 功能：返回A的每个元素的反正割值。,11双曲余割csch函数与反双曲余割asch函数 格式：A=csch(X) 功能：对X的每个元素求双曲余割。 格式：Xacsch(A) 功能：返回A的每个元素的反双曲余割值。,12余切cot函数与反余切acot函数 格式：A=cot(X) 功能：对X的每个元素求余切。 格式：Xacot(A) 功能：返回A的每个元素的反余切值。,1
29、3双曲余切coth函数与反双曲余切函数acoth函数 格式：A=coth(X) 功能：对X的每个元素求双曲余切。 格式：Xacoth(A) 功能：返回A的每个元素的反双曲余切值,指数、对数、幂运算,1指数exp函数与自然对数log函数 格式：Aexp(X) 功能：返回X每个元素的以e为底的指数值。 格式：Xlog(A) 功能：返回A的每个元素的自然对数。,2常用对数logl0函数 格式：Xlogl0(A) 功能：对A的每个元素求常用对数。,3 log2函数和幂函数pow2 格式：Xlog2(A) 功能：对A的每个元素计算其以2为底的对数。 格式：A=pow2(X) 功能：A的元素为2对X每个元
30、素求幂得到。,4平方根函数 格式：A=sqrt(X) 功能：返回对数组X每个元素求平方根的数组，如果元素为负或复数，则结果为复数。,复数的基本运算,1求复数实部函数real 格式：X=real(Z) 功能：返回复数Z的每个元素的实部。 2求复数的虚部函数imag 格式：Y=imag(Z) 功能：返回复数Z的每个元素的虚部。 3绝对值和复数模函数abs 格式：A=abs(Z) 功能：返回Z每个元素的绝对值。如果Z的元素是复数，则返回其模。,4求相角函数angle 格式：P=angle(Z) 功能：对Z的每个元素求相角。 5共轭函数conj 格式：C=conj(Z) 功能：对Z的每个元素求共轭复数
31、。,数据的取舍与保留,1向0方向取整数函数fix 格式：I=fix(X) 返回X中每个元素的最靠近零的整数。若X为复数，则分别对X的实部和虚部取整。,2向-方向取整数函数floor 格式：I=floor(X) 返回X中每个元素的最靠近该元素的最小整数。若X为复数，则分别对X的实部和虚部取最小整。,3向-方向取整数函数ceil 格式：I=ceil(X) 返回X中每个元素的最靠近该元素的最大整数。若X为复数，则分别对X的实部和虚部取最大整数。 4四舍五入取整数函数round 格式：I=round(X) 返回X中每个元素的最靠近该元素的整数。若X为复数，则分别对X的实部和虚部取整数。,5模数余函数mod 格式：M=mod(X，Y) 返回X关于Y的余数。6除后余数函数rem 格式：M=rem(X，Y) 功能：返回X关于Y的余数。,7符号函数sign 格式：s=sign(X) 如果X大于零，则返回1；如果X等于零，则返回0；如果X小于零，则返回-1。 如果X为复数，则算法为：SIGN(X) = X ./ ABS(X)。,
