第2讲matlab的数值分析

1第二讲MATLAB的数值分析2-1矩阵运算与数组运算矩阵运算和数组运算是MATLAB数值运算的两大类型，矩阵运算是按矩阵的运算规则进行的，而数组运算则是按数组元素逐一进行的。因此，在进行某些运算（如乘、除）时，矩阵运算和数组运算有着较大的差别。在MATLAB中，可以对矩阵进行数组运算，这时是把矩阵视为数组，运算按数组的运算规则。也可以对数组进行矩阵运算，这时是把数组视为矩阵，运算按矩阵的运算规则进行。1、矩阵加减与数组加减矩阵加减与数组加减运算效果一致，运算符也相同，可分为两种情况：（1）若参与运算的两矩阵（数组）的维数相同，则加减运算的结果是将两矩阵的对应元素进行加减，如A=[111;222;333];B=A;A+Bans=222444666（2）若参与运算的两矩阵之一为标量（1*1的矩阵），则加减运算的结果是将矩阵（数组）的每一元素与该标量逐一相加减，如A=[111;222;333];A+2ans=3334445552、矩阵乘与数组乘（1）矩阵乘矩阵乘与数组乘有着较大差别，运算结果也完全不同。矩阵乘的运算符为“*”，运算是按矩阵的乘法规则进行，即参与乘运算的两矩阵的内维必须相同。设A、B为参与乘运算的两矩阵，C为A和B的矩阵乘的结果，则它们必须满足关系Cm×n=Am×kBk×n。因此，参与运算的两矩阵的顺序不能任意调换，因为A*B和B*A计算结果很可能是完全不一样的。如：A=[111;222;333];B=A;2A*Bans=666121212181818F=ones(1,3);G=ones(3,1);F*Gans3G*Fans=111111111（2）数组乘数组乘的运算符为“.*”，运算符中的点号不能遗漏，也不能随意加空格符。参加数组乘运算的两数组的大小必须相等（即同维数组）。数组乘的结果是将两同维数组（矩阵）的对应元素逐一相乘，因此，A.*B和B.*A的计算结果是完全相同的，如：A=[11111;22222;33333];B=A;A.*Bans=111114444499999B.*Aans=111114444499999由于矩阵运算和数组运算的差异，能进行数组乘运算的两矩阵，不一定能进行矩阵乘运算。如A=ones(1,3);B=A;A.*Bans=111A*A???Errorusing==3Innermatrixdimensionsmustagree.3、矩阵除与数组除矩阵除分为矩阵右除和矩阵左除两种情况。矩阵右除的运算符为“/”，设A、B为两矩阵，则“A/B”是指方程X*B=A的解矩阵X。显然，矩阵右除运算对参与运算两矩阵的维数是有一定要求的，即矩阵A和B的列数必须相等。如A=[1111;2222;3333];B=[1111];X=A/BX=12X*Bans=11112222矩阵右除允许参与右除运算的矩阵B为标量，这时矩阵右除运算的结果是将矩阵A的每一元素逐一与该标量进行除法运算。如：A=[2468；8642]；B=2；A/Bans=12344321矩阵左除运算符为“\”，设A、B为两矩阵，则“A\B”是指方程B*X=A的解矩阵X。显然，矩阵左除运算对参与运算两矩阵的维数也有一定要求的，即矩阵A和B的行数必须相等。如数组右除的运算符为“./”，左除的运算符为“.\”。数组右除和左除的运算结果是完全等效的。设A、B为两同维矩阵，则“A./B”的运算结果是将矩阵A的每一个元素与矩阵B的对应元素相除。应注意的是，参与数组运算的两矩阵（数组）的大小必须相等。A=[223344;112233;445566];B=[123322;111111;225533];A./Bans=221122112233221122B./Aans=0.50000.50001.00001.00000.50000.500041.00001.00000.50000.50000.33330.33330.50000.50001.00001.00000.50000.50004、常用的矩阵运算函数（1）用size()函数计算矩阵A的维数，调用格式：d=size(A)%将矩阵A的行数和列数赋给变量d[m,n]=size(A)%将矩阵A的行数赋给变量m、列数赋给变量n（2）用rand()函数产生随机矩阵，，调用格式：rand(n)%产生值在0~1之间随机分布的n*n的随机方阵rand(m,n)%%生值在0~1之间随机分布的n*m的随机矩阵（3）计算矩阵长度（列数）的函数length()，调用格式：a=length(B)%将矩阵B的列数赋值给变量a（4）矩阵元素的求积运算函数prod()，调用格式：prod(A)%若A为向量，将计算矩阵A所有元素之积；若A为矩阵，将产生一行向量，其元素分别为矩阵A的各列元素之积。prod(A,k)%将对矩阵A按k定义的方向进行示积运算，若k=1则按列的方向求积，若k=2则按行的方向求积。（5）矩阵元素的求和运算函数sum()，调用格式同prod()函数。2-2多项式及其运算MATLAB提供了标准多项式的常用函数，包括求根、相乘、相除等。1．多项式的表达与创建MATLAB采用将多项式系数按幂次序排列形式的行向量来表征一多项式。设多项式为：0111asasasasAnnnn)(则表征该多项式的行向量表示为：][011aaaaPnn。因此，在MATLAB中，创建多项式即可创用建行向量的方法，直接输入按顺序排序的多项式系数即可。如，输入语句：A=[1221];即表示创建多项式12223sss，并赋值给变量A。2．多项式求根函数roots()用于对多项式求根，调用格式：p=roots(A)其中A为表征多项式的行向量，p返回该多项式的根（用列向量表示）。如B=[132];%创建多项式232ss5roots(B)%求多项式的根ans=-2-1A=[1221];roots(A)ans=-1.0000-0.5000+0.8660i-0.5000-0.8660i3．由指定根求多项式函数poly用于由给定根求多项式系数向量，调用格式为：A=poly(p)其中p为多项式根（行或列向量表示），A为返回的多项式系数（行向量表示）。如：p=[21];poly(p)ans=1-32可见roots()与poly()互为逆运算。4．多项式相乘（卷积）函数conv()用于求两个多项式的乘积多项式，调用格式：R=conv(A,B)其中A、B分别为表征两个多项式的行向量，R为返回的每间积多项式的系数向量（按降幂次序排列）。如A=[132]；B=[121]；R=conv(A,B)R=159725．多项式相除（解卷）函数deconv()用于进行两个多项式的相除运算，它是相乘运算（conv）的逆运算，其调用格式为：[B,t]=deconv(R,A)其中R为被除多项式，A为除数多项式，B为商多项式，t为余多项式。即多项式R除以多项式A后得商多项式B和余多项式t。如：R=[15972];A=[132];[B,t]=deconv(R,A)6B=121t=00000若多项式系数向量为零向量，则表示R能被A除尽。2-3线性代数方程（组）求解我们习惯将上组方程式以矩阵方式表示如下：AX=B其中A为等式左边各方程式的系数项，X为欲求解的未知项，B代表等式右边之已知项，要解上述的联立方程式，我们可以利用矩阵左除\做运算，即是X=A\B。如果将原方程式改写成XA=B其中A为等式左边各方程式的系数项，X为欲求解的未知项，B代表等式右边之已知项。注意上式的X,B已改写成列向量，A其实是前一个方程式中A的转置矩阵。上式的X可以矩阵右除/求解，即是X=B/A。若以反矩阵运算求解AX=B,X=B，即是X=inv(A)*B，或是改写成XA=B,X=B，即是X=B*inv(A)。我们直接以下面的例子来说明这三个运算的用法：A=[32-1;-132;1-1-1];%将等式的左边系数键入B=[105-1]';%将等式右边之已知项键入，B要做转置X=A\B%先以左除运算求解X=%注意X为行向量-256C=A*X%验算解是否正确C=%C=B105-1A=A';%将A先做转置B=[105-1];X=B/A%以右除运算求解的结果亦同X=%注意X为列向量105-1X=B*inv(A);%也可以反矩阵运算求解72-4信号的MATLAB表示举例1、指数信号atAeMATLAB中表示：y=A*exp(a*t)%example2-1decayingexponentialsignalA=1;a=-0.4;t=0:0.001:10;yt=A*exp(a*t);plot(t,yt)2、正弦信号)sin(tA0MATLAB中表示：y=A*sin(w0*t+phi)%example2-2sinusoidalsignalA=1;w0=2*pi;phi=pi/6;t=0:0.001:8;yt=A*sin(w0*t+phi);plot(t,yt)3、抽样函数)(tSaMATLAB中抽样函数用sinc函数表示：y=sinc(t)%example2-3Samplefunctiont=-3*pi:pi/100:3*pi;yt=sinc(t/pi);plot(t,yt)4、矩形脉冲信号MATLAB中矩形脉冲信号用rectpuls函数表示：y=rectpuls(t,width)%width缺省值为1%example2-4产生一幅度为1，宽度为width以零点为对称的矩形波t=0:0.001:4;T=1;yt=rectpuls(t-2*T,T);plot(t,yt)8axis([0,4,-0.5,1.5])5、三角波脉冲信号MATLAB中三角波脉冲信号用tripuls函数表示：y=tripuls(t,width,skew)用以产生一幅度最大为1，宽度为width的三角波，非零范围（-width/2,width/2）。skew定义skew=2tmax/width，tmax为三角波顶点坐标，skew取值：-1~+1，决定三角波的形状，缺省时为0。%example2-5t=-4:0.001:4;yt=tripuls(t,4,0.5);plot(t,yt)6、周期矩形脉冲信号y=square(w0*t,duty_cycle);用以产生一个幅度为+1和-1，基波频率为w0的矩形脉冲信号。duty_cycle指一个周期内正脉冲的宽度和负脉冲宽度的百分比，缺省值为1。%example2-6squarewavet=0:0.0001:5;A=1;T=1;w0=2*pi/T;ft=A*square(w0*t,20);plot(t,ft)axis([0,5,-1.5,1.5])7、周期三角波信号y=sawtooth(w0*t,width);产生一个基波频率为w0的三角波信号，其幅度为+1和-1，-1出现在nT点。设在[0,T]区间内+1出现的位置为tmax，则width=tmax/T。width缺省值为0.5.%example2-7triangularwavet=0:0.0001:5;A=1;T=1;w0=2*pi/T;yt=sawtooth(w0*t,1);plot(t,yt)axis([0,5,-1.5,1.5])98、信号的尺度变换、翻转、平移将信号x(t)写成函数，函数名为x2_1，程序如下：functionyt=x2_1(t)yt=0.5*(t+2).*(t=-2&t=0)+1*(t0&t=1);调用函数x2_1，即可现出x(t)波形，实现x(-2t-3)：%example1-8triangularwavet=-3:0.001:2;subplot(2,1,1)plot(t,x2_1(t))title(‘x(t)’)axis([-3,2,-1,2])subplot(2,1,2)plot(t,x2_1(-2*t-3))title(‘x(-2t-3)’)axis([-3,2,-1,2])9、连续信号的微分与积分连续信号的微分可由diff近似计算：h=0.001;x=0:h:pi;y=diff(sin(x.^2))/