系统仿真工具MATLAB.ppt

系统仿真工具MATLAB, 在仿真数值求解中的应用 东北大学信息学院智能所 郝培锋,MATLAB基本内容,MATLAB语言简介 程序设计基础 MATLAB数值线性代数计算 数值微积分及方程求解 系统工具箱(略),一 MATLAB语言简介,1 MATLAB具有交互式的语言环境，同时可以进行正常的程序化设计。 2 最早的语言主要是为进行数值代数及矩阵计算设计的系统，随着系统计算和绘图功能的完善其在各领域得到极其广泛的应用。 3 与其他语言有相应的接口，编程简单，函数功能齐全。 4 各种高级算法、工具箱等，Simulink面向对象的程序设计仿真环境等。,二 MATLAB程序设计基础,1 直接点击MATLAB图标可以自动建立交互式语言环境，基本的赋值语句结构 赋值变量 = 赋值表达式 例：A=1,2,3;7,8,0 表示将一个矩阵赋给A 函数调用 返回变量列表=函数名(输入变量列表) 例：x,y=meshgrid(-3:0.1:3,-2:0.2:2) 生成x,y的一个平面分割，横坐标在-3与3，纵坐标在-2与2间步长都为0.1。,二 MATLAB程序设计基础,2 简单的矩阵运算 矩阵转置 B=A 加减法运算(同阶矩阵) C1=A+B C2=A-B 乘法运算(可相乘) C=A*B 矩阵左除 C=AB 为AX=B的解 矩阵右除 C=A/B 为XA=B的解 矩阵翻转 C=fliplr(A), C=flipud(A), C=rot90(A) 矩阵乘方 C=A3 矩阵点运算 C=A.*B C=A.B C=A./B,二 MATLAB程序设计基础,3 MATLAB基本编程方式和程序语句 循环结构 for 循环变量= s1:s2:s3 循环体语句组 end 例: mysum=0; for i=1:1:100, mysum=mysum+i; end; mysum 运行结果: mysum=5050 相应的循环结构: while 逻辑变量 循环体语句组 end,二 MATLAB程序设计基础,3 MATLAB基本编程方式和程序语句 条件转移结构 if 逻辑变量 条件块语句组 end 例: mysum=0; for m=1:1:1000, if(mysum10000),break; end mysum=mysum+m; end; m,mysum 运行结果: m mysum = 142 10011,二 MATLAB程序设计基础,4 MATLAB函数编写规范 函数基本结构 Function 返回变量表=函数名(输入变量表) % 注释说明语段 输入、返回变量格式检测 函数体语句 注意：函数中有两个默认的输入输出变量个数检测变量nargin和nargout。(详细例题请参阅相关指导书),三 MATLAB数值线性代数计算,1 特殊矩阵的输入与表示 零、幺和单位阵 A=zeros(m,n), B=ones(m,n), C=eye(m,n) 随机元素阵、对角矩阵 A=rand(m,n), B=diag(a1,an),或B=diag(V) Vandermonde矩阵 例： C=1,2,3,4,5; V=wander(C) 结果：,1 1 1 1 1 16 8 4 2 1 625 125 25 5 1,三 MATLAB数值线性代数计算,2 特殊矩阵特征参数运算 矩阵行列式 A=2,4,0;-3,1,-12;0,-1,1, B=det(A) 矩阵迹与秩 B=tr(A) C=rank(A) 矩阵范数 N=norm(A,选项n)， 选项n对于范数的定义 n=1 表示矩阵的1-范数，即A=max n=2 矩阵的2-范数，即A= 而 表示A的最大奇异值的平方。,三 MATLAB数值线性代数计算,3 矩阵的特征多项式、特征方程与特征根 特征多项式：构造sI-A并求其行列式 C(s)=det(sI-A)=sn+c1sn-1+cn-1s+cn 特征多项式函数 C=poly(A) 如 A=1 2 3;4 5 6;7 8 0; B=poly(A) B= 1.000 -6.000 -72.000 -27.000 特征值与特征向量：v,d=eig(A),v= 0.7471 0.2998 0.2763 -0.6582 0.7075 0.3884 0.0931 0.6400 0.8791,d=-0.3884 0 0 0 12.1229 0 0 0 -5.7345,三 MATLAB数值线性代数计算,4 矩阵的逆 非奇异矩阵求逆： B=inv(A) 矩阵广义逆： 任何矩阵A都有 B=Pinv(A) 5 矩阵的非线性运算 面向矩阵元素的运算对象是矩阵的每一个元素。常用的一些面向元素的非线性函数有： abs() asin() acos() aqrt() log() log10() exp() real() imag() conj() sin() cos() round() floor() ceil(),三 MATLAB数值线性代数计算,面向矩阵元素的运算 例：A=1,2,3;4,5,6;7,8,0; exp(A) ans = 1.0e+003* 0.0027 0.0074 0.0201 0.0546 0.1484 0.4034 1.0966 2.9810 0.0010 log(A) warning:log of zero. ans = ,三 MATLAB数值线性代数计算,面向整个矩阵的非线性运算 矩阵指数的运算 expm(A),用Taylor级数展开eA,对矩阵进行运算： eA = = I+A+1/2!A2+1/3!A3+1/m!Am+ MATLAB还允许对矩阵求对数logm()、平方根sqrtm()和任意函数funm()。其中funm()是通过直接编写函数实现任意函数的整个矩阵运算。 例：funm(A,sin) 得对A的sin值，而它是面向矩阵元素的函数。,三 MATLAB数值线性代数计算,还可以通过直接编写函数的方法求得A的sin函数值。 例：function E=sinm1(A) E= zeros(size(A); F=A; k=1; while norm(E+F-E,1)0 E=E+F; F=-A2*F/(k+1)*(k+2);k=k+2; end E=sinm1(A) 实验表明funm(A,sin)与E结果是完全一致的。,四 数值微积分及方程求解,1 数值差分与积分运算 数值差分计算 差分函数diff(),设向量yi,则处理后的新向量为yi+1-yi,显然新的比原向量维数少1。 例如:A=2,3,1;0,-10,5;-3,8,20;1,1,3; diff(A) ans = -2 -13 4 -3 18 15 -4 -7 -17,四 数值微积分及方程求解,数值积分算法 如果定积分函数为: 则积分调用格式： y,n=quad(F,a,b,tol) 可以通过inline()函数定义该积分函数，如 f=inline(1/sqrt(2*pi)*exp(-x.2/2),x); y,kk=quad(f,-8,8) y = 0.99999997742635 kk= 860 其中积分表示式为 通常用有限区间长代替无限积分区间,四 数值微积分及方程求解,通过编写函数实现数值积分 Function y=myquad(x) y=1/sqrt(2*pi)*exp(-x.2/2); 还是对上述的无限区间上的定积分有： format long;y,kk=quad(myquad,-8,8) y = 0.99999997742635 kk= 504 也可以增大积分区间长进行计算 y,kk=quad(myquad,-15,15) y = 1.00000007062479 kk=860,四 数值微积分及方程求解,多(双)重定积分数值求解 一般形式的表示式 基本调用格式：y=dblquad(函数名,xM,xm,yM,ym) 例如： 编写函数。 Function z=my2dfun(x,y global kk;kk=kk+1; z=exp(-x.2/2).*sin(x.2+y); Y=dblquad(my2dfun,-2,2,-1,1) Y=1.57456866245358 kk=1781,四 数值微积分及方程求解,2 常微分方程数值解法 四阶定步长Runge-Kutta法 记tk时刻状态向量xk,定义四个附加向量： K1=h*f(tk,xk) K2=h*f(tk+h/2,xk+k1/2) K3=h*f(tk+h/2,xk+K2/2) K4=h*f(tk+h/2,xk+K3) 因而tk+h时刻的四阶R-K数值解由下式得出： xk+1=xk+(K1+2K2+2K3+K4)/6,四 数值微积分及方程求解,Runge-Kutta法函数 Function tout,yout=rk4(odefile,tspan,y0) t0=tspan(1);th=tspan(end); if length(tspan)=3, h=tspan(3); else, h=tspan(2)-tspan(1);end tout=t0:h:th;yout=; for t=tout, k1=h*eval(odfil(t,y0); k4=h*eval(odfil(t+h,y0+k3); y0=y0+(k1+2*k2+2*k3+k4)/6;yout=yout;y0; end,四 数值微积分及方程求解,4/5阶变步长Runge-Kutta算法函数 对于4/5阶的算法定义了一个误差向量(略)，根据它变换计算步长，所以称为自适应变步长法。 Matlab提供了ode23()和ode45()两个采用变步长的R-K算法，调用格式为： t,x=ode23(方程函数名，t0,tf,x0,选项) t,x=ode45(方程函数名，t0,tf,x0,选项) 注意方程函数名