计算机仿真实验报告 MATLAB

注：以下所有程序均在MATLAB7.0下运行通过。

实验一MATLAB语言编程一、实验目的：熟悉MATLAB语言及其环境，掌握编程方法。

要求认真听取实验指导老师讲解与演示。

二、具体实验内容、步骤、要求：1.运行交互式学习软件，学习MATLAB语言2.在MATLAB的命令窗口下键入如下命令：INTRO(注意：intro为一个用MATLAB语言编写的幻灯片程序，主要演示常用的MATLAB语句运行的结果。

)然后，根据显示出来的幻灯片右图按钮进行操作，可按START→NEXT→NEXT按钮一步步运行，观察。

3.自编程序并完成上机编辑，调试，运行，存盘：（1）、用MATLAB命令完成矩阵的各种运算。

例如：求出下列运算结果，并上机验证。

（1）A(:,1) %取矩阵A的第一列元素ans =11213141（2）A(2,:) %取矩阵A的第二行元素ans =21 22 23 24（3）A(1:2,2:3) %取矩阵A第一二行第二三列的元素ans =12 1322 23（4）A(2:3,2:3) %取矩阵A第二三行第二三列的元素ans =22 2332 33（5）A(:,1:2) %取矩阵A第一列与第二列元素ans =11 1221 2231 3241 42（6）A(2:3) %取矩阵A第二行与第三行的首列元素ans =21 31（7）A(:) %将矩阵A的所有元素按一列排列ans =11213141122232421323334314243444（8）A(:,:) %显示矩阵Aans =11 12 13 1421 22 23 2431 32 33 3441 42 43 44（9）ones(2,2) %建立一个两行两列的全1矩阵ans =1 11 1（10）eye(2) %建立一个二维的单位矩阵ans =1 00 1（2）、绘制数学函数的图形：例如：，理解数组运算与矩阵运算的功能。

MATLAB程序如下：t=0:0.1:8; %建立向量ty=1-2*exp(-t.*sin(t)); %计算向量t的函数向量yplot(t,y); %利用plot命令绘图xlabel('t');ylabel('y=1-2*e^(-t*sin(t))'); %注释横坐标与纵坐标图1.1 所对应的函数图像4.理解函数文件与命令文件的区别，并自编函数文件并调用。

命令文件可以理解为简单的M文件，命令文件中的变量都是全局变量。

函数文件是在命令文件的基础之上多添加了一行函数定义行，其代码组织结构和调用方式与对应的命令文件截然不同。

函数文件是以函数声明行“function...”作为开始的，其实质就是用户往MATLAB函数库里边添加了子函数，函数文件中的变量都是局部变量，除非使用了特别声明。

函数运行完毕之后，其定义的变量将从工作区间中清除。

而命令文件只是将一系列相关的代码结合封装，没有输入参数和输出参数，即不自带参数，也不一定要返回结果。

而多数函数文件一般都有输入和输出变量，并见有返回结果。

举例如下：已知非线性方程组求系统在初始条件为的数值解。

1）定义方程组的m函数equ.mfunction xs=equ(xi)x=xi(1);y=xi(2);z=xi(3);xs=zeros(3,1);xs(1)=sin(x)+y^2+log(z)-7;xs(2)=3*x+2^y-z^3+1;xs(3)=x+y+z-5;2）方程组求解x0=[1 1 1];xyz=fsolve(@equ,x0)xyz =0.5991 2.3959 2.0050在调用fsolve命令求解方程时调用了自定义的equ函数。

5.学会通过Help，熟悉MATLAB中为用户提供的功能各异的函数与命令。

只需在MATLAB命令窗口中键入help 命令名或函数名，按回车键后即会显示该命令或函数的相关说明。

根据说明即可自学如何使用该命令或函数。

例如键入help eig实验二数值积分算法练习与函数调用一、实验目的：理解数值积分法，熟练掌握MATLAB的函数调用。

二、实验示例介绍：1.用Euler法求初值问题的数值解：设方程如下：取步长h=0.1，名为FZSYZ1.M上机用如下程序名FZSYZ1.M可求出数值解。

2.在MATLAB中提出了现成的数值积分函数，如ode1、ode23、ode45求解微分方程组，下面介绍ode23它为二阶/三阶的RKF方法在MATLAB的ToolBox。

文件夹中的MATLAB/funfun下的M文件，在此介绍其调用方法与应用示例如下：[t,x]=ode23(‘系统函数名’，)其中，系统函数名为描述系统状态方程的M函数的名称，该函数名在调用时应该用引号括起来（文件名字串），为起始与终止时间，为系统初始变量的值（列向量），为控制解的精度（缺省值为=在ode23中），trace为输出形式控制变量，非零则程序运行的每步都显示出来。

trace为缺省值时——不显示中间结果。

t为输出参数返回积分的时间离散值（列向量）。

X：输出参量，返回每个时间点值的解的列向量。

注意：系统函数的编写格式为固定的。

例如：若MATLAB求解初值问题的解，其方程如下：第一步：编写如下程序：function xdot=fun21(t,x)xdot=x-t^2并以fun21.m存盘。

第二步：编写如下程序并以fzsy22.m存盘。

t0=0;tf=3;tol=1e-6;x0=1;trace=1;[t,x]=ode23('fun21',t0,tf,tol,trace)plot(t,x)第三步：在命令窗口下运行fzsy21即可求出x的解，并画出曲线。

3.实验具体内容、步骤、要求：（1）运行交互式软件中函数调用，学习程序，（2）试将（2-2）方程改为用Euler编程求解，试比较用ode23求解结果。

MATLAB源程序如下：%ode23求解t0=0;tf=3;tol=1e-6;x0=1;trace=1;[t,x]=ode23('fun21',t0,tf,tol,trace);plot(t,x,'r');hold on;%Euler求解h1=0.1;t=[t0:h1:tf];n=length(t);u=x0;uu(1)=u;for i=2:ndu=u-t(i-1)^2;u=du*h1+u;uu(i)=u;endplot(t,uu,'b--');hold off;图2.1 式（2-2）用Euler编程与ode23算法求解结果（3）试将（2-1）方程改为用ode23算法调用求解，并试比较结果。

图2.2 式（2-1）用Euler编程与ode23算法求解结果从比较可知：Euler需要小步长，不能自动确定步长；ode23算法自动确定步长，如果误差小于指定的精度，进行下一步，否则缩小步长，因此ode23算法比Euler算法更为准确，但占用的计算时间更长。

（4）利用ode23或ode45求解线性时不变系统微分方程并绘制出的曲线。

式中MATLAB程序清单如下：%函数(y(t) ) ̇=Ay(t)function ys=fun22(t,y)ys=zeros(2,1);ys(1)=-0.5*y(1)+y(2);ys(2)=-y(1)-0.5*y(2);%利用ode45求解其数值解并作图t0=0;tf=4;tol=1e-6;y0=[0;1];trace=1;y=zeros(2,1);[t,y]=ode45('fun22',t0,tf,y0,tol,trace);plot(t,y(:,1),'r');hold on;plot(t,y(:,2),'b');hold off；图2.3 第(4)问利用ode45算法求解所得y(t)曲线（5）求出与的单位阶跃响应，并分别求出状态空间模型。

MATLAB程序清单如下：num1=[2];den1=[1 2 1];g1=tf(num1,den1);%构造连续传函G1(S)num2=[1];den2=[2 3 3 1];g2=tf(num2,den2);%构造连续传函G2(S)step(g1);hold on;%画出G1(S)的阶跃响应step(g2);hold off;%画出G2(S)的阶跃响应[A1,B1,C1,D1]=tf2ss(num1,den1)%得到G1(S)的状态空间描述[A2,B2,C2,D2]=tf2ss(num2,den2)%得到G2(S)的状态空间描述从MATLAB窗口中得到的结果可知，对应的状态空间模型为：对应的状态空间模型为：图2.4 传递函数G1(S)与G2(S)对应的单位阶跃响应选做题一：已知系统传递函数为求对应的零极点模型，绘制系统阶跃响应。

MATLAB程序清单如下：num=[200 400];den=conv([1 1],[1 10 42]);g=tf(num,den);%构造连续传递函数[z,p,k]=tf2zp(num,den)%将普通传函转化为零极点模型step(g);在MATLAB命令窗中显示如下结果：z =-2p =-5.0000 + 4.1231i-5.0000 - 4.1231i-1.0000k =200所以对应的零极点模型为。

图2.5 选做题一中的传递函数对应的单位响应曲线图选做题二：已知一个二阶系统为，分别取R=0与R=0.5时，系统模型取不同步长的数值解（h=0.5~0.001）。

MATLAB程序清单如下：function ys=cstep(t,y)R=0;ys=zeros(2,1);ys(1)=y(2);ys(2)=-2*R*y(2)-y(1);[t,y]=ode45(@cstep,[0:0.5:5],[100 0]);plot(t,y(:,1),'r--');hold on;tol=1e-6;trace=1;t0=0;tf=5;[t1,y1]=ode45('cstep',t0,tf,[100 0],tol,trace);plot(t1,y1(:,1),'b');hold off在MATLAB命令窗中显示如下结果：图2.6 R=0时的仿真图2.7 R=0.5时的仿真从图2.6与图2.7可以看出，步长越短，仿真结果越接近准确值。

当步长取0.5时，可以看到与步长取0.001时有明显偏离，而采取自动步长的ode45算法与步长取0.001的仿真结果一致。

（6）设方程式为用各种数值积分法与不同步长求方程式的数值解，并比较之。

MATLAB程序清单如下：%欧拉法求解子函数function [uu,t]=euler(t,h1,y0)n=length(t);u=y0;uu(1)=u;for i=2:ndu=-40*u;u=du*h1+u;uu(i)=u;end%主程序t0=0;tf=0.5;y0=2;h1=0.01;t1=[0:h1:1];[uu1,t1]=euler(t,h1,y0);%步长取0.01，用欧拉法求解plot(t1,uu1,'r');hold on;h2=0.001;%步长取0.01，用欧拉法求解t2=[0:h2:1];[uu2,t2]=euler(t2,h2,y0);plot(t2,uu2,'b');hold on;tol=1e-6;trace=1;u0=1;[t,u]=ode23('fun26',t0,tf,y0,tol,trace);%用ode23求解，步长自定plot(t,u,'g--');hold off;gridgtext('Euler0.01');gtext('Euler0.001');gtext('ode23');在MATLAB命令窗中显示如下结果：图2.8 采用不同算法对第（6）题微分方程求解结果示意图从图2.5可以看出，采用欧拉法求解时，步长选取越短，其结果越接近真实值，选用龙阶-库塔法时，由于程序自动变步长直到达到预设精度，因此利用该方法求得的结果较准确。