MATLAB仿真技术.ppt

1、第三章 Matlab在科学运算中的应用,系统仿真本身涉及了大量的数学运算，计算机有巨大的求解数学运算的潜力，其求解问题的能力主要取决于算法与算法的代码实现。Matlab是一种高效的、高精度的科学运算语言，用它可以轻易地求解看似复杂的数学问题，且较容易上手。,3.1 解析解与数值解,现代科学与工程的进展离不开数学。数学家们感兴趣的问题和其他科学家、工程技术人员所关注的问题是不同 的。数学家往往对数学问题的解析解，或称闭式解(closed-form solution)和解的存在性严格证明感兴趣，而工程技术人员一般对如何求出数学问题的解更关心，也就是，能用 某种方法获得问题的解则是工程技术人员更关心

2、的问题。而获得这样解的最直接方法就是通过数值解法技术。在实际应用中，至少有两种情况需要数值解法： 解析解不存在时 解析解存在但不实用时,3.2数值线性代数问题及求解 3.2.1矩阵的特征参数运算 1.矩阵的行列式(determinant) Matlab提供了内建函数det()，得用它可以直接求出矩阵的行列式。 A=1 2 3;4 5 6;7 8 0;det(A) ans= 27 2.矩阵的迹(trace) 假设一个方阵A=aij，则矩阵A的迹定义为该矩阵对角线上各个元素之和。由代数理论可知，矩阵的迹和矩阵的特 征值之和是相同的，矩阵的迹可由trace()函数求出。,3.矩阵的秩(rank) 若

3、矩阵所有的列向量中共有rc个线性无关，则称矩阵的列秩为rc ，如果rc=m，则称A为列满秩矩阵。相应地，若矩阵所有的行向量中共有rr个线性无关，则称矩阵的行秩为rr ，如果rr=n，则称A为行满秩矩阵。可以证明，矩阵的行秩和列秩是相等的，故称之为矩阵的秩，记作： rank(A)= rc= rr 矩阵的秩也表示矩阵中行列式不等于0的子式的最大阶次。所谓子式就是从原矩阵中任取k行k列所构成的子矩阵。 Matlab提供了一个内建函数rank(A,)来用数值方法求取一个已知矩阵A的数值秩，其中为机器精度。如无特殊说明可用rank(A)求出A矩阵的秩。,4.矩阵的特征值和特征向量 对一个矩阵A来说，如果

4、存在一个非零的向量x，且有一个标量满足 Ax= x 则称为A矩阵的一个特征值，而x为对应于特征值的特征向量。矩阵的特征值与特征向量由Matlab提供的函数eig()可以容易的求出，调用格式为V,D=eig(A)，其中A为要处理的矩阵，D为一个对角矩阵，其对角线上的元素为矩阵A的特征值，而每个特征值对应的V矩阵的列为该特征值的特征向量，该矩阵是一个满秩矩阵。Matlab的矩阵特征值的结果满足AV=VD，且每个特征向量元素的平方和均为1。如果调用该函数时只给出一个返回变量，则将只返回矩阵A的特征值。即使A为复数矩阵也可由该函数求出其特征值和特征向量。, A=1 2 3 ; 4 5 6;7 8 0;

5、v,d=eig(A) v = -0.29982462710385 -0.74706733429712 -0.27625410987281 -0.70747178057871 0.65820191579521 -0.38842553979250 -0.63999130671192 -0.09306253848733 0.87909570970133 d = 12.12289378463239 0 0 0 -0.38838384240732 0 0 0 -5.73450994222508 eig(A) ans = 12.12289378463239 -0.38838384240732 -5.734

6、50994222508,5.矩阵的特征多项式、特征方程和特征根 构造一个矩阵sI-A，并求出该矩阵的行列式，则可得 出一个多项式C(s) C(s)=det(sI-A)=sn+c1sn-1+cn-1s+cn 这样的多项式C(s)称为矩阵A的特征多项式，其中系数 ci, i=1,2,n称为矩阵的特征多项式系数。 Matlab提供了求取特征多项式的函数C=poly(A)，而 返回的C为一个行向量，其各个分量为矩阵A的降幂排列的 特征多项式系数。该函数的另外一个调用格式是：如果给 定A为向量，则假定该向量是一个矩阵的特征根，由此求出 该矩阵的特征多项式系数。,令特征多项式等于0所构成的方程称为该矩阵的

7、特征方程，而特征方程的根称为该矩阵的特征根。可以调用Matlab函数V=roots(C)而直接获得。其中V为特征方程式的解，即原矩阵的特征根。 A=1 2 3;4 5 6;7 8 0; B=poly(A) B = 1.0000 -6.0000 -72.0000 -27.0000 roots(B) ans = 12.1229 -5.7345 -0.3884,3.2.2.多项式及多项式矩阵的求值 1.多项式的加减乘除 一个n次的多项式可以表示成： p(x)=anxn+an-1xn-1+a1x+a0 因此在Matlab中，可以用一个长度为n+1的行向量来表示p(x)如下： p=an an-1 a1

8、a0 多项式的加减，可直接由向量的加减而推出。 P1+P2 P1-P2 注意：矩阵P1与P2的长度必须一致。,多项式的乘除，可使用conv及deconv命令来完成 P3=conv(p1,p2) q,r=deconv(p1,p2) 结果中q为商式，r为余式。,2.多项式求值 多项式的求值可由polyval()函数直接完成，对于多项式矩阵来说，则可以由polyvalm()函数来完成，这两个函数的调用格式为C=polyval(p,x)或B=polyvalm(p,A)，其中p为多项式系数降幂排列构成的向量，x为一个向量，而A为一个给定的矩阵，这时返回的矩阵B为下面的矩阵多项式的值 B=a1An+a2A

9、n-1+anA+an+1I 其中I为和A同阶次的单位矩阵。 P=1 2 1; X=0:0.5:3; Y=polyval(P,X) plot(X,Y,o),若要计算P(A)，A为一方阵，可用polyvalm命令如下： A=1 2；3 4 B=polyvalm(P,A) 此结果和B=A2+2*A+1是一样的。 若上式改为B=polyval(P,A)，则其结果和B=A.2+2*A+1一样。这一点需注意。 多项式的根可以用roots命令求得。,3.部分分式展开 命令residue可用于计算一个分式的部分分式展开。若A(s)和B(s)为多项式，且B(s)无重根，则分式A(s)/B(s)可表示为： 其中p

10、1,p2,pn为A(s)的根(或是A(s)/B(s)的极点)， r1,r2,rn为常数，K(s)为一多项式，例如求3s+8/s2+5s+6的部分分式展开，可输入如下： B=3 8;A=1 5 6; r,p,k=residue(B,A) Residue命令也可将部分分式转回原来的形式，如： B2,A2=residue(r,p,k),3.2.3 线性代数问题的解析求解 Matlab不但提供了丰富的线性代数问题的数值解求解函数，还和著名的符号运算语言Maple有机结合，并以其为符号运算引擎提供了符号运算工具箱，该工具箱的函数可以用来求解线性代数问题的解析解。 在使用符号运算工具箱之前，需要把一些变量

11、声明为“符号变量”，以区别于常规的数值变量。将x声明为符号变量的语句格式为： x=sym(x,变量其它说明) 其中，“变量其它说明”选项包括实数real、正数positive等，对一般变量没有必要声明这样的类型。矩阵变量可以由sym(A)语句声明。还可以使用syms命令同时声明若干个变量及其类型，如：,syms a b c d syms a b c d positive,3.2.4 微积分问题的解析解,limit() 极限 diff() 微分 int() 不定积分 taylor() Taylor幂级数展开 syms x, limit(3*x2-x+1)/(2*x2+x+1)(x3/(1-x),

12、syms x;y=1/(x2+4*x+3)*sin(x);pretty(y) y1=diff(y,x) simple(y1) pretty(ans) y2=int(y1,x) y2=simple(y2) pretty(y2) syms a b;y3=int(y1,x,a,b) simple(y3) pretty(ans) y4=taylor(y,x,20) syms n;symsum(2/2n+2/3n,n,1,inf) %无穷级数求解,例1：图示电路，R=5、Ra=25 、Rb=100 、Rc=125 、Rd=40 、Re=37.5 ,求图中40V直流电压源的输出电流。,1、建立模型,R1=

13、(Rb*Rc)/(Ra+Rb+Rc) R2=(Ra*Rb)/(Ra+Rb+Rc) R3=(Ra*Rc)/(Ra+Rb+Rc),2、根据模型编写M文件进行仿真（eg1.m）,clear; V=40;R=5;Ra=25;Rb=100;Rc=125;Rd=40;Re=37.5; R1=(Rb*Rc)/(Ra+Rb+Rc); R2=(Ra*Rb)/(Ra+Rb+Rc); R3=(Ra*Rc)/(Ra+Rb+Rc); Req=R+R1+1/(1/(R2+Re)+1/(R3+Rd); i=V/Req,3、根据模型编写网孔电流方程进行仿真(eg2.m),clear; syms V R Ra Rb Rc Rd

14、 Re; r11=R+Rb+Rd; r22=Ra+Rb+Rc; r33=Ra+Rd+Re; A=r11 -Rb -Rd;-Rb r22 -Ra;-Rd -Ra r33; B=40 0 0; I=inv(A)*B I=simple(I) I1=I(1) I1=subs(I1,V,R,Ra,Rb,Rc,Rd,Re,40,5,25,100,125,40,37.5),3.3 Matlab下的常微分方程求解函数 Matlab提供了两个常微分方程求解的函数ode23()和ode45()。这两个函数分别采用了二阶三级的RKF方法和四阶五级的RKF方法，并采用自适应变步长的求解方法，即当解的变化较慢时采用较大

15、的计算步长，从而使得计算速度快，当方程的解变化较快时，积分步长会自动地变小，从而使得计算的精度很高。这两个函数和调用格式为： t,x=ode23(方程函数名,tspan,x0,选项,附加参数) t,x=ode45(方程函数名,tspan,x0,选项,附加参数),其中，选项可以通过odeget()和odeset()函数来设 置，具体的常用选项如下： RelTol为相对误差容许上限，默认值为0.001(即0.1%的相对误差)，在一些特殊的微分方程求解中，为了保证较高的精度，还应适当减小该值。 AbsTol为一个向量，其分量表示每个状态变量允许的绝对误差，其默认值为10-6。当然也可以自由设置其值，

16、以改变求解精度。 MaxStep为求解方程最大允许步长 在一般应用中没有必要修改其默认属性。这里用到的 方程函数名为描述系统状态方程的M函数名称，应该用引 号括起来。,变量tspan一般为仿真范围，例如取tspan=t0,tf，其中t0和tf分别 为用户指定的起始和终止计算时间，从5.0版本开始允许totf。另外， tspan还可以取作不等间距的时间向量。变量x0是系统初始状态变量 的值，注意，应该使该向量的元素个数等于系统状态变量的个数，否 则将给出错误信息。如果想采用默认的零初始状态，则可以在该处给 出一个空矩阵。 有了这些参数，就可以调用这两个函数对系统直接进行求解了。 此函数返回两个变

17、量t和x，其中t为求解的时间变量，因为采用变步长 的求解算法，所以得出的t向量并不一定是等间隔的；另一个变量x为 状态变量在各个时刻所组成的列向量所构成矩阵的转置。有了t和x这 两个变量，在求解过程结束后即可以用plot(t,x)来绘制出解的结果曲 线。,这里要用到的方程函数名的编写格式是固定的，如果 其格式没有按照要求去编写，则将得出错误的求解结果。 方程函数的引导语句为： function xdot=方程函数名(t,x,flag,附加参数) 其中t为时间变量，x为方程的状态变量，而xdot为状 态变量的导数。注意，即使微分方程是非时变的，也应该 在函数输入变量列表中写上t占位。可见，如果想

18、编写这样 的函数，首先必须已知原系统的状态方程模型。 如果有附加参数需要传递，则可以将其在原函数中给 出，若有多个附加参数，则它们之间应用适号隔开，且应 该确保它们与主调函数完全对应。另外应该用一个变量flag 来占位。,3.3.1常微分方程举例 通过几个例子演示Matlab中常微分方程求解的方法， 并指出求解过程中可能遇到的问题和这些问题的解决方 案。 Lorenz模型的状态方程为： 若令其初值x1(0)=x2(0)=0,x3(0)=10-10，则,考虑著名的Van der Pol方程 使用“函数句柄”的概念，其函数编写的格式与前面介绍的一致，所不同的是，实际应用中若需要带附加参数，则在写函

19、数文件时不用flag变量占位，这样在ode45()时，则不用引用函数名，而直接用句柄即可。,例：假设给定的隐式微分方程如下,如果能证明A(x)为非奇异的矩阵，则直接就能将该方程变换为标准的一阶微分方程组的形式，即x=A-1(x)B(x)，套用各种Matlab函数求解对应的方程。事实上，由于Matlab能较好的处理奇异问题，所以可以信赖Matlab来求解矩阵的逆，看是否有奇异的错误信息提示，如果没有相应的错误信息则应该相信得出的结果。,3.4 刚性方程的Matlab求解 在许多领域中，经常遇到一类特殊的常微分方程，其中一些解变化缓慢，另一些变化快，且相差较悬殊，这类方程常常称为刚性方程，又称Stiff方程