第2讲1 MATLAB的数值分析.docx

第2讲1 MATLAB的数值分析.docx

《第2讲1 MATLAB的数值分析.docx》由会员分享，可在线阅读，更多相关《第2讲1 MATLAB的数值分析.docx（14页珍藏版）》请在冰豆网上搜索。

第2讲1MATLAB的数值分析

第2讲

（1）MATLAB的数值分析

2-1矩阵运算与数组运算

矩阵运算和数组运算是MATLAB数值运算的两大类型，矩阵运算是按矩阵的运算规则进行的，而数组运算则是按数组元素逐一进行的。

因此，在进行某些运算（如乘、除）时，矩阵运算和数组运算有着较大的差别。

在MATLAB中，可以对矩阵进行数组运算，这时是把矩阵视为数组，运算按数组的运算规则。

也可以对数组进行矩阵运算，这时是把数组视为矩阵，运算按矩阵的运算规则进行。

1、矩阵加减与数组加减

矩阵加减与数组加减运算效果一致，运算符也相同，可分为两种情况：

（1）若参与运算的两矩阵（数组）的维数相同，则加减运算的结果是将两矩阵的对应元素进行加减，如

A=[111;222;333];

B=A;

A+B

ans=

222

444

666

（2）若参与运算的两矩阵之一为标量（1*1的矩阵），则加减运算的结果是将矩阵（数组）的每一元素与该标量逐一相加减，如

A=[111;222;333];

A+2

ans=

333

444

555

2、矩阵乘与数组乘

（1）矩阵乘

矩阵乘与数组乘有着较大差别，运算结果也完全不同。

矩阵乘的运算符为“*”，运算是按矩阵的乘法规则进行，即参与乘运算的两矩阵的内维必须相同。

设A、B为参与乘运算的两矩阵，C为A和B的矩阵乘的结果，则它们必须满足关系Cm×n=Am×kBk×n。

因此，参与运算的两矩阵的顺序不能任意调换，因为A*B和B*A计算结果很可能是完全不一样的。

如：

A=[111;222;333];

B=A;

A*B

ans=

666

121212

181818

F=ones（1,3）;

G=ones（3,1）;

F*G

ans

3

G*F

ans=

111

111

111

（2）数组乘

数组乘的运算符为“.*”，运算符中的点号不能遗漏，也不能随意加空格符。

参加数组乘运算的两数组的大小必须相等（即同维数组）。

数组乘的结果是将两同维数组（矩阵）的对应元素逐一相乘，因此，A.*B和B.*A的计算结果是完全相同的，如：

A=[11111;22222;33333];

B=A;

A.*B

ans=

11111

44444

99999

B.*A

ans=

11111

44444

99999

由于矩阵运算和数组运算的差异，能进行数组乘运算的两矩阵，不一定能进行矩阵乘运算。

如

A=ones（1,3）;

B=A;

A.*B

ans=

111

A*A

?

?

?

Errorusing==>

Innermatrixdimensionsmustagree.

3、矩阵除与数组除

矩阵除分为矩阵右除和矩阵左除两种情况。

矩阵右除的运算符为“/”，设A、B为两矩阵，则“A/B”是指方程X*B=A的解矩阵X。

显然，矩阵右除运算对参与运算两矩阵的维数是有一定要求的，即矩阵A和B的列数必须相等。

如

A=[1111;2222;3333];

B=[1111];

X=A/B

X=

1

2

X*B

ans=

1111

2222

矩阵右除允许参与右除运算的矩阵B为标量，这时矩阵右除运算的结果是将矩阵A的每一元素逐一与该标量进行除法运算。

如：

A=[2468；8642]；

B=2；

A/B

ans=

1234

4321

矩阵左除运算符为“\”，设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./B

ans=

221122

112233

221122

B./A

ans=

0.50000.50001.00001.00000.50000.5000

1.00001.00000.50000.50000.33330.3333

0.50000.50001.00001.00000.50000.5000

4、常用的矩阵运算函数

（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采用将多项式系数按幂次序排列形式的行向量来表征一多项式。

设多项式为：

则表征该多项式的行向量表示为：

。

因此，在MATLAB中，创建多项式即可创用建行向量的方法，直接输入按顺序排序的多项式系数即可。

如，输入语句：

A=[1221];

即表示创建多项式

，并赋值给变量A。

2．多项式求根

函数roots（）用于对多项式求根，调用格式：

p=roots（A）

其中A为表征多项式的行向量，p返回该多项式的根（用列向量表示）。

如

B=[132];%创建多项式

roots（B）%求多项式的根

ans=

-2

-1

A=[1221];

roots（A）

ans=

-1.0000

-0.5000+0.8660i

-0.5000-0.8660i

3．由指定根求多项式

函数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=

15972

5．多项式相除（解卷）

函数deconv（）用于进行两个多项式的相除运算，它是相乘运算（conv）的逆运算，其调用格式为：

[B,t]=deconv（R,A）

其中R为被除多项式，A为除数多项式，B为商多项式，t为余多项式。

即多项式R除以多项式A后得商多项式B和余多项式t。

如：

R=[15972];

A=[132];

[B,t]=deconv（R,A）

B=

121

t=

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为行向量

-2

5

6

C=A*X%验算解是否正确

C=%C=B

10

5

-1

A=A';%将A先做转置

B=[105-1];

X=B/A%以右除运算求解的结果亦同

X=%注意X为列向量

105-1

X=B*inv（A）;%也可以反矩阵运算求解

2-4信号的MATLAB表示举例

1、指数信号

MATLAB中表示：

y=A*exp（a*t）

%example2-1decayingexponentialsignal

A=1;a=-0.4;

t=0:

0.001:

10;

yt=A*exp（a*t）;

plot（t,yt）

2、正弦信号

MATLAB中表示：

y=A*sin（w0*t+phi）

%example2-2sinusoidalsignal

A=1;w0=2*pi;phi=pi/6;

t=0:

0.001:

8;

yt=A*sin（w0*t+phi）;

plot（t,yt）

3、抽样函数

MATLAB中抽样函数用sinc函数表示：

y=sinc（t）

%example2-3Samplefunction

t=-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）

axis（[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-5

t=-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-6squarewave

t=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-7triangularwave

t=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]）

8、信号的尺度变换、翻转、平移

将信号x（t）写成函数，函数名为x2_1，程序如下：

functionyt=x2_1（t）

yt=0.5*（t+2）.*（t>=-2&t<=0）+1*（t>0&t<=1）;

调用函数x2_1，即可现出x（t）波形，实现x（-2t-3）：

%example1-8triangularwave

t=-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））/h;

定积分可由quad函数或quadl函数实现:

quad（‘function_name’,a,b）%function_name为被积函数名，a和b指积分区间

例：

求example2-5所示三角波求微分和积分

先将三角波写成MATLAB函数，名为x2_2:

functionyt=x2_2（t）

yt=tripuls（t,4,0.5）;

用diff和quad函数可实现三角波的微分和积分

%example2-9differentiation

h=0.001;x=-3:

h:

3;

y1=diff（x2_2（t））*1/h;

plot（t（1:

length（t）-1）,y1）

title（'dx（t）/dt'）

%example1-10integration

t=-3:

0.1:

3;

forn=1:

length（t）

y2（n）=quad（'x2_2',-3,t（n））;

end

plot（t,y2）

title（'integralofx（t）'）

2-5练习

1、实现下列连续时间信号

（1）

（画t=-1~3）

（2）

（3）

（画t=-1~3）

（4）

（画t=-1~3）

提示：

画阶跃信号可先创建如下函数，并保存为Heaviside.m，再调用即可。

functionf=Heaviside（t）

f=（t>0）;

调用格式：

f=Heaviside（t）

2、已知x（t）的波形如下图，

（1）编写实现该信号波形的函数x（t）的MATLAB函数；

（2）画出x（t），x（0.5t）和x（2-0.5t）。

答案：

t=-4:

0.01:

10;

subplot（3,1,1）;

plot（t,x（t））;title（'x（t）'）;

axis（[-4,10,-2,3]）

t=-4:

0.01:

10;

subplot（3,1,2）;

plot（t,x（0.5*t））;title（'x（0.5t）'）;

axis（[-4,10,-2,3]）

t=-6:

0.01:

8;

subplot（3,1,3）;

plot（t,x（2-0.5*t））;title（'x（2-0.5t）'）;

axis（[-6,8,-2,3]）