2MATLAB的数值计算分析

第二讲MATLAB的数值计算——matlab具有出色的数值计算能力，占据世界上数值计算软件的主导地位数值运算的功能创建矩阵矩阵运算多项式运算线性方程组数值统计线性插值函数优化微分方程的数值解一、命令行的基本操作创建矩阵的方法直接输入法规则：矩阵元素必须用[]括住矩阵元素必须用逗号或空格分隔

在[]内矩阵的行与行之间必须用分号分隔

矩阵元素可以是任何matlab表达式，可以是实数，也可以是复数，复数可用特殊函数I，j输入

a=[123;456]x=[2pi/2;sqrt(3)3+5i]

矩阵元素符号的作用逗号和分号的作用

逗号和分号可作为指令间的分隔符，matlab允许多条语句在同一行出现。

分号如果出现在指令后，屏幕上将不显示结果。注意：只要是赋过值的变量，不管是否在屏幕上显示过，都存储在工作空间中，以后可随时显示或调用。变量名尽可能不要重复，否则会覆盖。当一个指令或矩阵太长时，可用•••续行冒号的作用

用于生成等间隔的向量，默认间隔为1。用于选出矩阵指定行、列及元素。循环语句2.用matlab函数创建矩阵空阵[]—matlab允许输入空阵，当一项操作无结果时，返回空阵。rand——

随机矩阵（rand（3）、rand（2，3））eye——

单位矩阵zeros——全部元素都为0的矩阵ones——全部元素都为1的矩阵3.特殊矩阵矩阵符号说明company酉矩阵gallery试验矩阵hadamard哈德马德矩阵hankel汉克尔矩阵hilb希尔伯特矩阵invhilb逆希尔伯特矩阵magic魔方矩阵pascal帕斯卡矩阵randn随机矩阵，元素服从正态分布rosser对称特征值实验矩阵vander范德蒙矩阵wilkinson威尔金森特征值测试矩阵

还有伴随矩阵、稀疏矩阵、魔方矩阵、对角矩阵、范德蒙等矩阵的创建，就不一一介绍了。注意：matlab严格区分大小写字母，因此a与A是两个不同的变量。

matlab函数名必须小写。4.矩阵的修改

直接修改可用键找到所要修改的矩阵，用键移动到要修改的矩阵元素上即可修改。指令修改可以用A(,)=来修改。例如a=[120;305;789]a=120305789a(3,3)=0a=120305780还可以用函数subs修改，matlab6.0还可用find函数修改实际应用中，常需要提取矩阵的某一部分，或者将多个矩阵合并起来——裁剪和拼接对矩阵的裁剪与拼接操作主要是通过冒号运算符（:）来实现（1）重新排列>>x=[123;456;789];>>x(3:-1:1,:)ans=7894561235.矩阵裁剪与拼接（2）提取行提取：>>x=[123;456;789];>>x(1,:)ans=123列提取：>>x=[123;456;789];>>x(:,1)ans=147部分提取：>>x=[123;456;789];>>x(1:2,2:3)ans=2356重复提取：>>x=[123;456;789];>>x(:,[1,1,1])ans=111444777x(:,[1;2;2])（3）删除>>x=[123;456;789];>>x(:,1)=[]x=235689>>size(x)ans=32

Matlab中只能删除矩阵的整行或整列，x(:,1)=[]将矩阵x的第一列删除，同样，x(1,:)=[]将矩阵x的第一行删除。（4）拼接行拼接：两个矩阵的行数必须相同>>x=[123;456;789];>>y=ones(3);>>z=[x,y]z=123111456111789111列拼接：两个矩阵的列数必须相同>>x=[123;456;789];>>y=ones(1,3);>>z=[x;y]z=1234567891116.矩阵函数7.矩阵分解矩阵加、减（＋,－）运算规则：相加、减的两矩阵必须有相同的行和列两矩阵对应元素相加减。允许参与运算的两矩阵之一是标量。标量与矩阵的所有元素分别进行加减操作。二、矩阵运算2.矩阵乘（）运算规则：A矩阵的列数必须等于B矩阵的行数标量可与任何矩阵相乘。a=[123;456;780];b=[1;2;3];c=a*bc=143223

d=[-1;0;2];f=pi*df=-3.141606.2832

矩阵除的运算在线性代数中没有，有矩阵逆的运算，在matlab中有两种矩阵除运算a^p——a自乘p次幂

方阵>1的整数3.矩阵乘方——a^p,p^a对于p的其它值,计算将涉及特征值和特征向量，如果p是矩阵，a是标量a^p使用特征值和特征向量自乘到p次幂；如a,p都是矩阵，a^p则无意义。a=[1,2,3;4,5,6;7,8,9];a^2ans=303642668196102126150※当一个方阵有复数特征值或负实特征值时，非整数幂是复数阵。a^0.5ans=

0.4498+0.7623i0.5526+0.2068i0.6555-0.3487i1.0185+0.0842i1.2515+0.0228i1.4844-0.0385i1.5873-0.5940i1.9503-0.1611i2.3134+0.2717iinv——

矩阵求逆det——

行列式的值eig——

矩阵的特征值diag——

对角矩阵

’

——

矩阵转置sqrt——

矩阵开方expm——矩阵的指数logm——矩阵的对数cond——矩阵的条件数rank——矩阵的秩norm——矩阵的范数4.矩阵的其它运算5.矩阵的一些特殊操作矩阵的变维

a=[1:12];b=reshape(a,3,4)c=zeros(3,4);c(:)=a(:)矩阵的变向

rot90:旋转;fliplr:上翻;flipud:下翻矩阵的抽取

diag:抽取主对角线;tril:抽取主下三角;

triu:抽取主上三角矩阵的扩展关系运算

关系符号意义<<=>>===~=小于小于或等于大于大于或等于等于不等于采用冒号运算符

a=初值:步长:终值；函数linspace创建

a=linspace(初值,终值,元素个数)函数logspace

a=logspace(初始指数,终结指数,元素个数)示例：

>>a=logspace(0,3,3)a=1.0e+003*0.00100.03161.00006.一维数组的特殊创建方法

数组运算指元素对元素的算术运算，与通常意义上的由符号表示的线性代数矩阵运算不同（1）数组加减(.+,.-)a.+ba.-b7.矩阵的数组运算对应元素相加减（与矩阵加减等效）（2）数组乘除(

，./，.\）ab——a，b两数组必须有相同的行和列两数组相应元素相乘。a=[123;456;789];b=[246;135;7910];a.*bans=281841530497290a=[123;456;789];b=[246;135;7910];a*bans=253746558510985133172a./b=b.\aa.\b=b./aa./b=b.\a—

都是b的元素被a的对应元素除a.\b=b./a—

都是a的元素被b的对应元素除例:a=[123];b=[456];c1=a.\b;c2=b./ac1=4.00002.50002.0000c2=4.00002.50002.0000——

给出a,b对应元素间的商.（4）数组乘方(.^)—

元素对元素的幂例:a=[123];b=[456];z=a.^2z=1.004.009.00z=a.^bz=1.0032.00729.00（5）向量函数

一般标量函数都可以用于数组，此时函数作用于数组的每个元素：

>>x=1:5;>>sin(x)ans=0.84150.90930.1411-0.7568-0.9589

但有些函数只有作用于向量时才有意义，故称其为向量函数：

max——最大值

min——最小值

sum——和

length——长度（size）

mean——平均值

mediam——中间值

prod——乘积

sort——从小到大排列向量函数也可以作用于矩阵，此时运算结果为一个行向量，行向量的每个元素是函数作用于矩阵的相应列向量所得到的结果（1）字符串Matlab中的字符串一般是ASCII值的数值数组，它作为字符串表达式进行显示一个字符串是由单引号括起来的文本，是由字符组成的数组

>>s='Iamateacher';>>ss=Iamateacher>>u=abs(s)u=73329710932973211610197991041011148.文本操作>>str1=s(5:10)str1=atea>>y=setstr(u)%函数setstry=Iamateacheru=s(10:-1:5)输出内容？u=s’输出内容？字符串I’mateacher如何输入？可以像数组一样连接字符串可以用disp函数打印字符串变量表示的字符串，如disp(x);字符串也可以有多个行，但每行必须也有相同数目的列数，如果长度不等，则要以空格补齐示例：把字符串由小写变为大写

>>s='iamateacher';>>fori=1:14,s(i)=setstr(s(i)-('a'-'A'));end,ss=IAMATEACHER如果改为：

>>s='iamateacher';>>fori=1:14,s(i)=setstr(s(i)-(‘a’-‘A’)),end,s请考虑输出如果改为：

>>s=‘Iamateacher';>>fori=1:14,s(i)=setstr(s(i)-(‘a’-‘A’));end,s请考虑输出（2）字符串转换函数意义abs字符串到ASCII转换dec2hex十进制到十六进制转换fprintf把格式化的文本写到文件中或显示屏上fscanf从文件中读入有格式数据hex2dec十六进制字符串转换成十进制数hex2num十六进制字符串转换成IEEE浮点数int2str整数转换成字符串lower字符串转换成小写num2str数字转换成字符串setstrASCII转换成字符串sprintf把数据格式化，写给一个字符串sscanf用格式控制，把字符串转换成数字str2mat字符串转换成一个文本矩阵str2num字符串转换成数字upper字符串转换成大写>>rad=2.5;>>t=[‘Acircleofradius’,rad]%是否正确fscanf调用格式：

A=fscanf(fid,format)[A,count]=fscanf(fid,format,size)fprintf调用格式：

fprintf(obj,'cmd')fprintf(obj,‘format’,‘cmd’)%obj可以省略sscanf调用格式：

A=sscanf(s,format)A=sscanf(s,format,size)sprintf调用格式：

[s,errmsg]=sprintf(format,A,...)控制符符号说明示例\b退格符sprintf(‘WhoamI?\b’)（WhoamI）\n回车换行符sprintf(‘WhoamI\n?’)\t横向走纸符

sprintf('WhoamI?\tadf')\\、’’、%%输出特殊字符sprintf(‘WhoamI?\\')f浮点输出符sprintf('%.2f',pi)3.14sprintf('%.0f',pi)3e指数输出符sprintf('%.2e',pi)3.14e+000sprintf('%.0e',pi)3e+000g定长输出符sprintf('%10.2g',pi)________3.1sprintf('%10.0g',pi)_________3（3）字符串函数函数意义blanks(n)返回一个n个零或空格的字符串deblank去掉字符串中后拖的空格findstr从一个字符串内找出字符串isletter字母存在时返回真值（返回值为矩阵）isspace空格字符存在时返回真值（返回矩阵）isstr是一个字符串，返回真值lasterr返回上一个所产生Matlab错误的字符串strcmp两个字符串相同，则返回真strrep用一个字符串代替另一个字符串strtok找出字符串中的第一部分s='Thisisagoodexample.';[token,rem]=strtok(s)token=Thisrem=isagoodexample.matlab语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。

f(x)=anxn+an-1xn-1+……+a0

可用行向量p=[anan-1

……a1+a0]表示poly——

产生特征多项式系数向量特征多项式一定是n+1维的特征多项式第一个元素一定是1三、多项式运算例:a=[123;456;780];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^3-6x^2-72x-272.roots——

求多项式的根a=[123;456;780];p=poly(a)p=1.00-6.00-72.00-27.00r=roots(p)r=12.12-5.73——显然r是矩阵a的特征值

-0.39当然我们可用poly令其返回多项式形式p2=poly(r)p2=1.00-6.00-72.00-27.00matlab规定多项式系数向量用行向量表示，一组根用列向量表示3.conv,convs多项式乘运算例:a(x)=x2+2x+3;b(x)=4x2+5x+6;c=(x2+2x+3)(4x2+5x+6)a=[123];b=[456];c=conv(a,b)=conv([123],[456])c=4.0013.0028.0027.0018.00p=poly2str(c,'x')p=4x^4+13x^3+28x^2+27x+184.deconv多项式除运算a=[123];c=[4.0013.0028.0027.0018.00]d=deconv(c,a)d=4.005.006.00[d,r]=deconv(c,a)余数c除以a后的整数5.多项式微分matlab提供了polyder函数多项式的微分。命令格式：polyder(p):求p的微分polyder(a,b):求多项式a,b乘积的微分例：a=[12345];poly2str(a,'x')ans=x^4+2x^3+3x^2+4x+5b=polyder(a)b=4664poly2str(b,'x')ans=4x^3+6x^2+6x+4四、代数方程组求解matlab中有两种除运算左除和右除。对于方程ax+b，a为an×m矩阵，有三种情况：当n=m时，此方程成为“恰定”方程当n>m时，此方程成为“超定”方程当n<m时，此方程成为“欠定”方程

matlab定义的除运算可以很方便地解上述三种方程1.恰定方程组的解方程ax=b(a为非奇异)x=a-1

b

矩阵逆两种解:x=inv(a)

b—

采用求逆运算解方程x=a\b—

采用左除运算解方程方程ax=ba=[12;23];b=[8;13];

x=inv(a)*b

x=a\bx=x=22.0033.00

=

ax=b例:x1+2x2=82x1+3x2=13左除、右除的异同2.超定方程组的解方程ax=b,m<n时此时不存在唯一解。方程解(a'a)x=a'bx=(a'

a)-1a'b——

求逆法

x=a\b——matlab用最小二乘法找一个准确地基本解。例:x1+2x2=12x1+3x2=23x1+4x2=3a=[12;23;34];b=[1;2;3];

解1x=a\b

解2x=inv(a'

a)

a'

b

x=x=1.001.0000.00

=

ax=b3.欠定方程组的解

当方程数少于未知量个数时,即不定情况,有无穷多个解存在。matlab可求出两个解：用除法求的解x是具有最多零元素的解是具有最小长度或范数的解，这个解是基于伪逆pinv求得的。x1+2x2+3x3=12x1+3x2+4x3=2a=[123;234];b=[1;2];x=a\bx=pinv(a)

bx=x=1.000.8300.330-0.17=ax=b五、微分方程求解微分方程求解的仿真算法有多种，常用的有Euler（欧拉法）、RungeKutta(龙格-库塔法。Euler法称一步法，用于一阶微分方程当给定仿真步长时：所以

yn+1=yn+h·f(xn,yn)n=0,1,2…y(x0)=y0RungeKutta法龙格-库塔法：实际上取两点斜率的平均斜率来计算的，其精度高于欧拉算法。龙格-库塔法：ode23ode45

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文件functionxdot=wf(t,x)xd