实验一面向微分方程的数值积分法仿真.docx

实验一 面向微分方程的数值积分法仿真
一、实验目的
1．掌握数值积分法的基本概念、原理及应用；
2．用龙格-库塔法解算微分方程，增加编写仿真程序的能力；
3．分析数值积分算法的计算步长与计算精度、速度、稳定性的关系；
4. 对数值算法中的“病态问题”进行研究。
二、实验内容
1、已知系统微分方程及初值条件

取步长，试分别用欧拉方程法和RK4法求时的值，并将求得的值与解析解比较(将三个解绘于同一坐标中，且用数值进行比较)，说明造成差异的原因。（①编程完成；②选用MATLAB ode函数完成。）
程序代码如下：
t0=0;
tf=2;
h=0.1;
y1=1;
y2=1;
y3=1;
t1=0;
t2=0;
t3=0
n=round(tf-t0)/h;
for i=1:n
y1(i+1)=y1(i)+h*(2*h+y1(i));
t1=[t1,t1(i)+h];
end
for i=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=[t2t2(i)+h];
end
for i=1:n
y3(i+1)=2*exp(t3(i))-t3(i)-1;
t3=[t3,t3(i)+h];
end
plot(t1,y1,'r',t2,y2,'g',t3,y3,'k')
实验结果如下；

分析：红线为用欧拉法得到的结果，绿线为用四阶龙格—库塔法得到的结果，蓝线为根据解析方程得到的结果。其差异原因主要有两个：1、二者的方法不同，欧拉法是根据一阶微分方程计算得到的，龙格—库塔法是根据四阶微分方程得到的；2、由于步长取为0.1，所以得到的图像与解析解之间存在差异，若将步长取小，则得到的解将更靠近解析解。
2、已知系统的传递函数为
在单位阶跃输入下，系统响应的解析解为

试分别用欧拉方程法和RK4法对系统进行仿真（编程完成）：
1）比较两种数值积分解与解析解得逼近程度；（绘图）
程序代码如下：
num=[40.6];den=[1 10 27 22.06];
[a,b,c,d]=tf2ss(num,den);
a=[0 1 0; 0 0 1; -22.06 -27 -10];
b=[0;0;1];
c=[40.6 0 0];

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;
for i=1:N
k1=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];
end
for i=1:N
x=X2(:,i)+h*(a*X2(:,i)+b*u);
y=c*x;
X2=[X2,x];
Y2=[Y2,y];
t2=[t2,t2(i)+h];
end
for i=1:N
y=1.84-4.95*i*exp(-1.88*i)-1.5*exp(-1.88*i)-0.34*exp(-6.24*i);
Y3=[Y3,y];
t3=[t3,t3(i)+h];
end
plot(t1,Y1,'r',t2,Y2,'g',t3,Y3,'b')

当h=0.01时的结果

当h=0.01时的结果

分析：这是我得到的结果，发现两个方法得到的结果与实际结果都有较大差距，当是龙格—库塔法更接近实际的结果。
2）改变步长，分析步长对数值解精度的影响；
改变步长后，发现只是两根仿真得到的曲线靠近了，但是与实际曲线仍然是差距很大，这是经过仔细的检查和讨论我觉得程序还是对的，不知道错在哪里了。
3）不断加大步长，分析计算稳定性的变化。
当取h=0.5时，得到的结果：

加大步长后结果得到的结果不稳定，不能够很好的对系统进行仿真，另外，由于系统步长选择偏大，根据解析解得到的结果也与实际值有了一定的差距，但是如果步长取得不一样又无法比较。
3、求下图所示系统的阶跃响应的数值解。（，，，，）分析、对系统响应的影响。（①编程用RK4求解；②Simulink）

程序代码如下：
k=1;
a=conv([1 0 0],conv([0.25 1],[0.25 1]))
b=[2*k k]
X0=[0 0 0 0];
V=1; n=4;T0=0;Tf=10; h=0.25;R=1;
b=b/a(1);a=a/a(1);A=a(2: n+1);
A=[rot90(rot90(eye(n-1,n)));-fliplr(A)];
B=[zeros(1,n-1),1]';
m1=length(b);
C=[fliplr(b),zeros(1,n-m1)];
Ab=A-B*C*V;
X=X0';y=0;t=T0;
N=round(Tf-T0)/h;
for i=1:N

k1=Ab*X+B*R;
k2=Ab*(X+h*k1/2)+B*R;
k3=Ab*(X+h*k2/2)+B*R;
k4=Ab*(X+h*k3)+B*R;
X=X+h*(k1+2*k2+2*k3+k4)/6;
y=[y,C*X];
t=[t,t(i)+h];
end
[t',y']
plot(t,y)
得到结果如下：
K=0.5,v=1时的结果如下：

K=1，v=1时，

K=1，V=0.5时

K=1,v=5时，

K=2，v=1时

分析：当k取值增大，v值不变时，系统输出的波头增多，而且也变陡，稳态精度降低，当k增加到一定程度时系统便发散了（即不稳定了）。当v值增大，k值不变时，波头也是变多变陡，当v值增大到一定程度时系统便不稳定了。
②Simulink

4、已知系统状态状态方程为

采用RK4法，步长分别取，求系统的零输入响应，并绘图分析各状态变量的响应状态及产生的原因。（提示：病态系统）
程序代码如下：
a=[-21 19 -20;19 -21 20;40 -40 -40];
x=[1;0;-1];X=x;t=0;
t0=t;tf=2;h=0.01;
n=round(tf-t0)/h;
for i=1:n
x=X(:,i)+h*a*X(:,i);
X=[X,x];
t=[t,t(i)+h];
end
b=X(1,:);c=X(2,:);d=X(3,:);
plot(t,b,'r',t,c,'g',t,d,'b')

当h=0.01时得到的结果

当h=0.02时得到的结果

当h=0.04时得到的结果

分析：如图，当h=0.01时，在t=0.2s以后系统输出便趋于平稳，当取h=0.02时，系统输出振荡剧烈，趋于稳定的时间也变长，当取h=0.04后，系统输出呈发散振荡形式。当h=0.06后系统仍然是发散的，即当h的取值改变时，原先稳定的系统变得不稳定了，这便是病态系统。但是这个结果与书上的不同。
三、 实验报告要求
记录完成实验内容所采取的步骤、方法和结果。回答思考题。
四、 预习要求及思考题
要求实验前，必须预习实验知识点，按实验内容的要求的确定仿真方案，完成程序设计，以便在实验时进行调试、分析。
思考题：1、在进行仿真计算时，是否选用的数值积分法的阶次越高越好？
答：阶次不是越高越好。（1）阶次越高计算公式也越复杂，每一步需要计算的次数也更多，需要的时间也更长。（2）每一阶的龙格库塔法都有对应的稳定区域，当hλ(λ为系统的特征方程根)超过稳定区域时，即使阶次高得到的结果也未必是稳定的。当hλ得值接近稳定边界时，误差也会增大。
2、选用数值积分法进行仿真的原则。
答：（1）精度 仿真结果的精度主要受三项误差的影响：
1）截断误差：由算法本身的精度阶次所决定。
2）舍入误差：由计算机字长决定。
3）累积误差：由以上两项误差随计算时间长短累积情况决定。
（2）计算速度 计算速度取决于所用的数值方法和计算步长。
（3）稳定性 数值稳定性主要与计算步长h有关，不同的数值方法对h都有不同的稳定性限制范围，且与被仿真对象的时间常数也有关系。
实验二
 第一部分  面向结构图的数值积分法仿真
一、实验目的
加深理解连续系统面向结构图仿真的原理及特点，进一步掌握数字积分法解算微分方程的方法，增加编写仿真程序的能力。

二、实验内容
1、用面向方框图的数字仿真方法对下列系统进行仿真。

2、求解下图所示系统在f=-1(t)阶跃扰动作用下第④、第⑤环节的动态过程。
分别用面向框图的数值积分法（RK4法）、MATLAB中有关系统建模的命令和Simulink三种方法求解。

第二部分  面向结构图的离散相似法仿真
一、实验目的
   1．掌握离散相似法的基本概念和原理；
2．典型环节离散系数的求取及差分方程表示；
3. 掌握非线性系统的数字仿真方法。
二、实验内容
1、已知控制系统结构图如图所示，设输入阶跃函数幅值Y0=10，滞环非线性参数s=1(滞环宽度），请用离散相似法编程和Simulink法对系统进行如下分析：
1）不考虑非线性环节影响时，求解y(t)的阶跃响应；
2）考虑非线性环节影响，其余参数不变，求解y(t)并与线性情况所得结果进行比较；
3）改变的滞环非线性参数s，分析该非线性对系统的影响。

P=[1 10 5 10;1 0.5 1 0;1 0.1 1 0;0 1 1 0];
WIJ=[1 0 1;2 1 1;3 2 1;4 3 1;1 4 -1];X0=[0 0 0 0];
Z=[0 0 0 0];S=[0 0 0 0];
h=0.01;L1=25;n=4;T0=0;Tf=10;nout=4;Y0=10;
function Uc=satu(Ur,S1)
if(abs(Ur)>=S1)
if(Ur>0) Uc=S1;
else Uc=-S1;
end
else Uc=Ur;
end
function Uc=dead(Ur,S1)
if(abs(Ur)>=S1)
if(Ur>0) Uc=Ur-S1;
else Uc=Ur+S1;
end
else Uc=0;
end
function [Uc,Ubb]=(Urb,Ur,Ucb,S1)
if(Ur>Urb)
if((Ur-S1)>=Ucb) Uc=Ur-S1;
else Uc=Ucb;
end
else if(Urlse Uc=Ucb;
end
else Uc=Ucb;
end
end
Ubb=Ur;
Y=zeros(n,1);X=Y;y=0;Uk=zeros(n,1);Ubb=Uk;
T=T0:h*L1:Tf;N=length(t);
for k=1:N-1
for l=1:L1
Ub=Uk;
Uk=W*Y+W0*Y0;
for i=1:n
if(Z(i)~=0)
if(Z(i)==1)
Uk(i)=satu(Uk(i),S(i));
end
if(Z(i)==2)
Uk(i)=dead(Uk(i),S(i));
if(Z(i)==2)
[Uk(i),Ubb(i)]=backlash(Ubb(i),Uk(i),Ub(i),S(i));
end
end
end
Udot=(Uk-Ub)/h;
Uf=2*Uk-Ub;

X=FI’.*X+FIM’.*Uk+FIJ’.*Udot;
Yb=Y;Y=FIC’.*X+FID’.*Uf;
for i=1:n
if(Z(i)~=0)
if(Z(i)==4)
Y(i)=satu(Y(i),S(i));
end
if(Z(i)==5)
Y(i)=dead(Uk(i),S(i));
if(Z(i)==6)
[Y(i),Ubb(i)]=backlash(Ubb(i),Y(i),Yb(i),Si));
end
end
end
end
y=[y,Y(nout)];
end

2、系统结构图如图所示，先理论分析该系统是否会产生自激振荡，若会求出振荡的振幅和频率，并用离散相似法编程和Simulink法验证分析的结果。

实验报告要求
     编程完成实验任务；记录完成实验内容所采取的步骤、方法；分析仿真结果。
实验三 采样控制系统的数字仿真
一、实验目的
1、掌握采样控制系统数字仿真的一般方法；
2、掌握采样周期与计算步长的关系与确定方法。
二、实验内容

1已知采样系统结构如图所示，分别用程序设计方法和Simulink法求系统的输出响应。（取仿真时间:Tf=10；采样周期:T=0.2；计算步长：h=0.01）(注：第一个框图中，分子分母中z的指数均为-1；第二个框图中e的指数为-Ts)

2、设某数字控制系统如图所示，采样周期为T=0.1s，初始状态。(注：第二个框图中e的指数为-Ts)

（1）数字控制器为PI调节器 ，其中
（2）数字控制器为PID调节器，其中

要求：
1）用程序设计法、MATLAB控制工具箱时域响应分析函数、Simulink法求系统在单位阶跃输入信号作用下的输出响应；
2）比较两种控制器作用下系统的响应，由此得出微分调节的作用。

三、完成实验报告

实验四 数字仿真技术的综合应用
一、实验目的
使学生在系统地学习了计算机仿真的基本概念和方法后，对计算机仿真的全过程有一个清楚的概念和认识，能够将所学的知识应用到实际系统中去。

二、实验内容

选用经典控制理论或现代控制理论中某一种你所熟知的控制器设计方法，为下列某一系统（三选一）进行控制器设计：
（一）一阶倒立摆系统
（二）水箱液位控制系统
（三）直流电动机调速系统

三、实验报告

内容要求完整， 包括建模、模型验证及分析、控制器设计、仿真实验等。其中仿真实验中要对系统所能达到的性能指标给予详细分析说明。