自控系统数字仿真实验

《自控系统数字仿真实验》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、实验一  面向微分方程的数值积分法仿真一、实验目的1．掌握数值积分法的基本概念、原理及应用；2．用龙格-库塔法解算微分方程，增加编写仿真程序的能力；3．分析数值积分算法的计算步长与计算精度、速度、稳定性的关系；4. 对数值算法中的“病态问题”进行研究。二、实验内容1、已知系统微分方程及初值条件取步长h=0.1 ，试分别用欧拉方程法和RK4法求t=2h时的y值，并将求得的值与解析解比较(将三个解绘于同一坐标中，且用数值进行比较)，说明造成差异的原因。（①编程完成；②选用MATLABode函数完成。）程序如下：t0=0;tf=2;h=

2、0.1;y1=1;y2=1;y3=1;t1=0;t2=0;t3=0;n=round(tf-t0)/h;fori=1:n%欧拉法y1(i+1)=y1(i)+h*(2*h+y1(i));t1=[t1,t1(i)+h];endfori=1:n%四阶龙格-库塔法k1=y2(i)+t2(i);k2=y2(i)+h*k1/2+t2(i)+h/2;k3=y2(i)+h*k2/2+t2(i)+h/2;k4=y2(i)+h*k3+t2(i)+h;y2(i+1)=y2(i)+h*(k1+2*k2+2*k3+k4)/6;t2=[t2,t2(i)+h];

3、endfori=1:ny3(i+1)=2*exp(t3(i))-t3(i)-1;t3=[t3,t3(i)+h];endplot(t1,y1,'r',t2,y2,'g',t3,y3,'k')分析:红线为用欧拉法得到的结果，绿线为用四阶龙格—库塔法得到的结果，蓝线为根据解析方程得到的结果。其差异原因主要有两个：（1）二者的方法不同，欧拉法是根据一阶微分方程计算得到的，龙格—库塔法是根据四阶微分方程得到的；（2）由于步长取为0.1，所以得到的图像与解析解之间存在差异，若将步长取小，则得到的解将更靠近解析解。 2、已知系统的传递函数为在单

4、位阶跃输入下，系统响应的解析解为 试分别用欧拉方程法和RK4法对系统进行仿真（编程完成）：1）比较两种数值积分解与解析解的逼近程度；（绘图）2）改变步长，分析步长对数值解精度的影响；3）不断加大步长，分析计算稳定性的变化。程序如下：a=[010;001;-22.06-27-10];b=[0;0;1];c=[40.600];X1=[0;0;0];t=0;Y1=0;X=0;u=1;Y2=0;Y3=0;X2=[0;0;0];x=0;h=0.1;t0=0;tf=2;t1=0;t2=0;t3=0;N=(tf-t0)/h;fori=1:Nk1

5、=a*X1+b;k2=b+a*(h*k1/2+X1);k3=b+a*(h*k2/2+X1);k4=b+a*(h*k3+X1);X1=X1+h*(k1+2*k2+2*k3+k4)/6;Y1=[Y1,c*X1];t1=[t1,t1(i)+h];endfori=1:Nx=X2(:,i)+h*(a*X2(:,i)+b*u);y=c*x;X2=[X2,x];Y2=[Y2,y];t2=[t2,t2(i)+h];endfori=1:Ny=1.84-4.95*i*exp(-1.88*i)-1.5*exp(-1.88*i)-0.34*exp(-6.

6、24*i);Y3=[Y3,y];t3=[t3,t3(i)+h];endplot(t1,Y1,'r',t2,Y2,'g',t3,Y3,'b')①当h=0.01时，仿真图如下：可见当步长很小时，数值积分解逼近解析解，其数值解得的精度也很高。②当h=0.1时，仿真图如下：可见随着步长的增加，其数值积分解与解析解的逼近程度减小，精度降低。③当h=0.3时，仿真图如下：可见其数值积分解已变得不稳定了，即步长不能取得太大3、求下图所示系统的阶跃响应y(t)的数值解。（v=1,k=1,开始时间t0=0,结束时间tf=10,h=0.25）分析k、

7、v对系统响应的影响。（①编程用RK4求解；②Simulink）(修正：G(s)表达式的分母中，0.25s+1改为（0.25s+1）的二次方)  程序如下：k=1;a=conv([100],conv([0.251],[0.251]))b=[2*kk];X0=[0000];v=1;n0=4;tf=10;h0=0.25;r=1;V=v;n=n0;T0=0;Tf=tf;h=h0;R=r;b=b/a(1);a=a/a(1);A=a(2:n+1);%首一化处理A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];%形

8、成能控标准型A阵B=[zeros(1,n-1),1]';%形成输入阵B(n维列向量)m1=length(b);%分子系数向量维数m+1C=[fliplr(b),zeros(1,n-m1)];%形成输出阵C(n维行向量)Ab=A-B*C*V;%形成闭