控制系统计算机仿真及辅助设计实验报告

控制系统计算机仿真及辅助设计实验报告实验一1-2用MATLAB绘制下列函数图像1.f(x){sin(x2)0≤x2ecos(x)2≤x4−lon(x)4≤x6解：编写程序x=0:0.01:6y=sin(x.^2).*(x=0&x2)+exp(cos(x)).*(x=2&x4)-log(x).*(x=4&x6)plot(x,y)图像2.f（x）{x2−y2|x−y|2e−|x−y||x−y|≫2编写程序x=linspace(-10,10)y=linspace(-10,10)[x,y]=meshgrid(x,y)z=(x.^2-y.^2).*(abs(x-y)2)+exp(-abs(x-y)).*(abs(x-y)=2)mesh(x,y,z)图像3.编写程序wn=sqrt(0.2)s=0.25/(wn*2)sys=2.5*tf([0.2],[10.250.2])t=linspace(0,100)impulse(sys,t)step(sys,t)单位脉冲响应图像单位阶跃响应图像实验二2-2.用MATLAB语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：1.G（s）=324327242410355024sssssss2..X=2.25-5-1.25-0.542.25-4.25-1.25-0.2520.25-0.5-1.25-121.25-1.75-0.25-0.750XuY=[0202]X1.解：（1）状态方程模型参数:编写MATLAB程序如下num=[172424];den=[110355024];[ABCD]=tf2ss(num,den)得到结果：A=-10-35-50-24100001000010,B=1000,C=172424,D=[0]所以模型为:.X=-10-35-50-24100001000010X+1000u,y=172424X(2)零极点增益:编写程序num=[172424];den=[110355024];[ZPK]=tf2zp(num,den)得到结果Z=-2.7306+2.8531,-2.7306-2.8531i,-1.5388P=-4,-3,-2,-1K=1(3)部分分式形式:编写程序num=[172424];den=[110355024];[RPH]=residue(num,den)得到结果R=4.0000,-6.0000,2.0000,1.0000P=-4.0000,-3.0000,-2.0000,-1.0000H=[]G(s)=46214321ssss阶跃响应程序:num=[172424];den=[110355024]sys=tf(num,den);t=linspace(0,10);step(sys,t);脉冲响应程序:num=[1,7,24,24];den=[1,10,35,50,24];sys=tf(num,den);t=linspace(0,10);impulse(sys,t);2.解：（1）传递函数模型参数：编写程序A=[2.25-5-1.25-0.52.25-4.25-1.25-0.250.25-0.5-1.25-11.25-1.75-0.25-0.75]；B=[4220]';C=[0202];D=[0];[numden]=ss2tf(A,B,C,D)得到结果:num=04.000014.000022.000015.0000den=1.00004.00006.25005.25002.2500324324s+14s+22s+15()s+4s+6.25s+5.25s+2.25Gs(2)零极点增益模型参数:编写程序A=[2.25-5-1.25-0.52.25-4.25-1.25-0.250.25-0.5-1.25-11.25-1.75-0.25-0.75]；B=[4220]';C=[0202];D=[0];[Z,P,K]=ss2zp(A,B,C,D)得到结果:Z=-1.0000+1.2247i-1.0000-1.2247i-1.5000P=-0.5000+0.8660i-0.5000-0.8660i-1.5000-1.5000K=4.0000表达式4s+1-1.2247is+1+1.2247i()s+0.5-0.866is+0.5+0.866is+1.5Gs（3）部分分式形式的模型参数：编写程序A=[2.25-5-1.25-0.52.25-4.25-1.25-0.250.25-0.5-1.25-11.25-1.75-0.25-0.75]；B=[4220]';C=[0202];D=[0];[numden]=ss2tf(A,B,C,D)[R,P,H]=residue(num,den)得到结果R=4.0000-0.00000.0000-2.3094i0.0000+2.3094iP=-1.5000-1.5000-0.5000+0.8660i-0.5000-0.8660iH=[]42.30942.3094()1.50.50.8660.50.866iiGsssisi脉冲响应程序:A=[2.25,-5,-1.25,-0.5;2.25,-4.25,-1.25,-0.25;0.25,-0.5,-1.25,-1;1.25,-1.75,-0.25,-0.75];B=[4;2;2;0];C=[0,2,0,2];D=0;Sys=ss(A,B,C,D);t=linspace(0,10);impulse(sys,t);阶跃响应程序A=[2.25,-5,-1.25,-0.5;2.25,-4.25,-1.25,-0.25;0.25,-0.5,-1.25,-1;1.25,-1.75,-0.25,-0.75];B=[4;2;2;0];C=[0,2,0,2];D=0;Sys=ss(A,B,C,D);t=linspace(0,10);step(sys,t);2-3.用欧拉法求下面系统的输出响应y(t)在0≤t≤1上，h=0.1时的数值。',(0)1yyy要求保留4位小数，并将结果与真解()tyte比较。解：欧拉法1'00*(,)(,)()kkkkkkyyhftyyftyyty（前向欧拉法,可以自启动）其几何意义：把f(t,y)在[,kkty]区间内的曲边面积用矩形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示（1）m文件程序为h=0.1;disp('函数的数值解为');%显示‘’中间的文字%disp('y=');%同上%y=1;fort=0:h:1m=y;disp(y);%显示y的当前值%y=m-m*h;end保存文件q2.m在matalb命令行中键入q2得到结果函数的数值解为y=10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487（2）另建一个m文件求解tye在t[0,1]的数值（%tye是',(0)1yyy的真解%）程序为h=0.1;disp('函数的离散时刻解为');disp('y=');fort=0:h:1y=exp(-t);disp(y);end保存文件q3.m在matalb命令行中键入q3函数的离散时刻解为y=10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679比较欧拉方法求解与真值的差别欧拉10.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679误差0-0.0048-0.0007–0.0118–0.0142–0.0160–0.0174–0.0183–0.0188-0.0192-0.0192显然误差与2h为同阶无穷小，欧拉法具有一阶计算精度，精度较低，但算法简单。2-4用二阶龙格库塔法求解2-3的数值解，并于欧拉法求得的结果比较。解：我们经常用到预报-校正法的二阶龙-格库塔法，112121'()2(,)(,)(,)kkkkkkhyykkkftykfthyhkftyy此方法可以自启动，具有二阶计算精度，几何意义：把f(t,y)在[,kkty]区间内的曲边面积用上下底为kf和1kf、高为h的梯形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示（1）m文件程序为h=0.1;disp('函数的数值解为');disp('y=');y=1;fort=0:h:1disp(y);k1=-y;k2=-(y+k1*h);y=y+(k1+k2)*h/2;end保存文件q4.m在matlab的命令行中键入q4显示结果为函数的数值解为y=10.90500.81900.74120.67080.60710.54940.49720.45000.40720.3685比较欧拉法与二阶龙格-库塔法求解.（误差为绝对值）真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679龙库10.90500.81900.74120.67080.60710.54940.49720.45000.40720.3685误差00.00020.00030.00040.00050.00060.00060.00060.00070.00060.0006明显误差为3h得同阶无穷小，具有二阶计算精度，而欧拉法具有以阶计算精度，二阶龙格-库塔法比欧拉法计算精度高。2-5．用四阶龙格-库塔法求解题2-3数值解，并与前两题结果相比较。解：四阶龙格-库塔法表达式112341213243(22)6(,)(,)22(,)22(,)kkkkkkkkkkhyykkkkkftyhhkftykhhkftykkfthyhk，其截断误差为5h同阶无穷小，当h步距取得较小时，误差是很小的.(1)编辑m文件程序h=0.1;disp('四阶龙格-库塔方法求解函数数值解为');disp('y=');y=1;fort=0:h:1disp(y);k1=-y;k2=-(y+k1*h/2);k3=-(y+k2*h/2);k4=-(y+k3*h);y=y+(k1+2*k2+2*k3+k4)*h/6;end保存文件q5.m在matlab命令行里键入q5得到结果四阶龙格-库塔方法求解函数数值解为y=10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679(2)比较这几种方法：对于四阶龙格-库塔方法真值10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679龙库10.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679误差00000000000显然四阶龙格-库塔法求解精度很高，基本接近真值。三种方法比较可以得到精度（四阶）〉精度（二阶）〉精度（欧拉）实验三3.2实验程序脉冲num=[0.8,0,-20];den=[1,0,-40,0];sys=tf(num,den);t=0:0.01:1;impulse(sys,t)实验图形阶跃num=[0.8,0,-20];den=[1,0,-40,0];