第3-1章-连续系统数值积分法仿真Matlab编程课件

第三章Matlab编程(数值积分法仿真)3.1连续系统数值积分法仿真编程思路目的：针对下面的系统编写通用的数值积分法仿真程序(3.1-1)对这样的系统进行仿真，实际上涉及到控制的计算、状态的数值积分计算和输出的计算。3个函数g,f,h确定后，就可以完整地描述一个系统。给定初值：u0，y0,系统中的向量都采用列向量表示1第三章Matlab编程(数值积分法仿真)3.1连续系统数function[t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)%函数功能：用数值积分法仿真%输入参数：tstart,tstop,h分别是起始时间、结束时间和仿真步长，是标量%x0,u0是状态、输入的初值，都是列向量%cnty是输出变量的个数%InteMethod时数值积分方法，可以是'EUL','RK2','RK4'%StateModel是状态模型的文件名%ControlFile是控制作用的文件名%OutputFile是系统输出的文件名%输出参数：t是仿真结果的时间序列%y是仿真结果系统的输出序列(1)数值积分法仿真框架函数2function[t,y]=w_DigiInteSimu(function[t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)t=[tstart:h:tstop];%t数一个行序列cntt=size(t,2);%返回列数y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果y0=eval([OutputFile,'(tstart,x0,u0)']);%计算初始输出y(:,1)=y0’;%将cury作为输出的第1列curx=x0;%当前一步的xcuru=u0;%当前一步的ucury=y0;%当前一步的yfori=1:1:cntt-1curu=eval([ControlFile,'(t(i),h,curx,curu,cury)']);%计算控制时传递的参数：当前时间，步长，当前状态和输出curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%单步积分运算cury=eval([OutputFile,'(t(i),curx,curu)']);%计算输出y(:,i+1)=cury‘;%将输出加入到输出序列里end3function[t,y]=w_DigiInteSimu(functionNewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel)%函数功能：单步积分运算,与模型方程有关%输入参数：curt是当前时间，h是数值积分步长% curx,curu分别是当前的状态和控制向量%InteMethod是积分算法，字符串类型，可以取'EUL','RK2','RK4'%StateModel是状态模型文件名称，字符串类型%输出参数：NewX是这一步计算的新的状态向量(2)单步数值积分法函数单步数值积分函数只是对微分方程组StateModel进行一步的计算，计算法由InteMethod参数指定，可以上欧拉法，RK2或RK4。4functionNewX=w_StepIntegral(cfunctionNewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel);ifInteMethod=='RK4'k1=eval([StateModel,'(curt,curx,curu)']);k2=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k1,curu)']);k3=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k2,curu)']);k4=eval([StateModel,'(curt+h,curx+h*k3,curu)']);NewX=curx+h*(k1+2*k2+2*k3+k4)/6;elseifInteMethod=='RK2'k1=eval([StateModel,'(curt,curx,curu)']);k2=eval([StateModel,'(curt+h,curx+h*k1,curu)']);NewX=curx+0.5*h*(k1+k2);else%欧拉法EULQk=eval([StateModel,'(curt,curx,curu)']);%eval用来执行一个函数，传递函数名和函数的输入参数NewX=curx+h*Qk;end5functionNewX=w_StepIntegral(c函数w_DigiInteSimu和w_StepIntegral构造了一个数值积分法仿真的框架，并不涉及具体的系统。具体的系统由StateModel,ControlFile,OutputFile参个参数决定，实际上就是三个函数文件名，这三个函数输入输出参数必须遵循特定的格式，在准备好由这3个函数描述的系统后，调用w_DigiInteSimu即可进行仿真。还需要一个调用w_DigiInteSimu进行仿真的程序，它指定模型文件，指定初始参数，并且对仿真结果绘图。6函数w_DigiInteSimu和w_StepIntegra3.2算法稳定性分析仿真编程针对下面的系统，求用解析法、欧拉法和RK4分别求解，计算欧拉法最大允许步长，将步长从0.1逐渐增大，比较三种解的效果。解：1）该系统是稳定的，解析解为2）用欧拉法计算本例时，其步长应该满足3）RK4步长条件式是一个高阶不等式，无法直接求解，只能用试探法确定RK4的步长上限。(3.2-1)73.2算法稳定性分析仿真编程针对下面的系统，求用解析法、欧我们使用3.1介绍的两个基本函数，同时要针对系统模型编写3个函数来描述该系统，最后编写一个实现仿真初值设置、仿真求解和仿真结果显示的函数。functionderX=w_SysStateEqu(curt,curx,curu)%functionderX=w_stateEqu(curt,curx,curu)%函数功能：在此函数里写系统的状态方程%输入参数：curt当前时间，curx和curu是当前状态和控制%输出参数：derX是状态的X的导数值derX(1)=-4*curx(1);（1）状态模型函数w_SysStateEqu在该函数里描述(3.2-1)式的状态微分方程8我们使用3.1介绍的两个基本函数，同时要针对系统模型编写3个functionOutputY=w_SysOutput(curt,curx,curu)%functionOutputY=w_SysOutput(curt,curx,curu)%函数功能：计算系统的输出向量%输入参数：curt是当前时间，curx,curu分别是当前的状态和控制向量%输出参数：OutputY是计算出来的输出向量OutputY(1)=curx(1);%根据具体的模型写输出方程（2）系统输出函数w_SysOutputfunctionControlU=w_SysControl(curt,h,curx,curu,cury)%ControlU=w_sysControl(curt,h,curx,curu,cury)%函数功能：计算系统的控制，如u=PID(t,curx,cury)%输入参数：curt是当前时间，h是数值积分步长％ curx,curu,cury分别是当前的状态、控制和输出向量%输出参数：ControlU是计算出来的控制向量ControlU=curu;%根据实际系统写控制（3）系统控制计算函数w_SysControl9functionOutputY=w_SysOutput(c（4）仿真组织函数functionsimu_Stability%函数功能：稳定性分析仿真tstart=0;tstop=7;h=0.7;StateModel='w_SysStateEqu';ControlFile='w_SysControl';OutputFile='w_SysOutput';x0=[1];u0=[0];cnty=1;[t,y1]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'EUL',StateModel,OutputFile,ControlFile);%用欧拉法求解[t,y4]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'RK4',StateModel,OutputFile,ControlFile);%用RK4求解y5=exp(-4*t);%计算解析解plot(t,y1,'--m',t,y2,'-.r',t,y5,'-k');%绘制叠加曲线xlabel('时间'); ylabel('y');title('算法稳定性分析仿真');该函数的主要任务是：1)设置仿真起始时间、结束时间、仿真步长、初始的x和u2)设置模型的输出方程、状态方程和控制方程所在函数文件名3)调用w_DigiInteSimu进行仿真计算4)将计算结果绘图进行比较。10（4）仿真组织函数functionsimu_Stabili修改simu_Stability函数里的参数h，就可以作出不同步长时的仿真结果。下面是步长从0.1不断增大时的仿真结果h=0.1时的仿真曲线。1】RK4法的曲线与解析解重合，说明RK4法非常准确。2】欧拉法有较大误差。11修改simu_Stability函数里的参数h，就可以作出不h=0.4时的仿真曲线。1】RK4法还是比较准确。2】欧拉法误差很大，出现衰减震荡。3】两种数值算法的步长条件仍然满足，算法仍然稳定。12h=0.4时的仿真曲线。12h=0.6时的仿真曲线。欧拉法的步长条件已经不满足，仿真解发散h=0.6时的仿真曲线。RK4法的步长条件仍然满足。但误差增大13h=0.6时的仿真曲线。欧拉法的步长条件已经不满足，仿真解发h=0.7时的仿真曲线。RK4法的步长条件已经不满足，仿真结果发散。14h=0.7时的仿真曲线。143.3卫星发射仿真卫星发射运动方程G为重力系数系统是高阶微分方程，不能直接用于仿真计算，需要转化为状态方程。定义状态变量：运动方程用极坐标表示，同时还需要用直角坐标输出。1.模型转换153.3卫星发射仿真卫星发射运动方程G为重力系数系统是高阶微得到系统仿真模型取X-Y平面坐标作为输出状态初值为：v是初始发射速度，可分别取该系统没有控制，主要是研究初始发射速度对轨道的影响。16得到系统仿真模型取X-Y平面坐标作为输出状态初值为：v是初始2.模型描述函数构造functionderX=sat_StateEqu(curt,curx,curu)G=401408;derX(1)=curx(3);derX(2)=curx(4);derX(3)=-G/(curx(1)*curx(1))+curx(1)*curx(4)*curx(4);derX(4)=-2*curx(3)*curx(4)/curx(1);derX=derX’;%转换为列向量(1)状态方程functionOutputY=sat_Output(curt,curx,curu)OutputY(1)=curx(1)*cos(curx(2));OutputY(2)=curx(1)*sin(curx(2));OutputY=OutputY’;%转换为列向量(2)输出方程(3)控制方程由于该系统没有控制，因而采用与3.2节同样的控制函数。172.模型描述函数构造functionderX=sat_S3.仿真组织函数functionsimu_Satelite%函数功能：对卫星发射轨道仿真tstart=0;StateModel='sat_StateEqu';ControlFile='w_SysControl';OutputFile='sat_Output'; u0=[0];cnty=2;h=200;tstop=100*h; v=10;x0=[6400,0,0,v/6400]‘;[t,y10]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解h=200;tstop=50*h; v=9;x0=[6400,0,0,v/6400]’;[t,y9]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解h=100;tstop=60*h; v=8;x0=[6400,0,0,v/6400]‘;[t,y8]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解plot(y10(1,:),y10(2,:),'-',y9(1,:),y9(2,:),'--',y8(1,:),y8(2,:),':');%绘制叠加曲线xlabel('x'); ylabel('y'); title('卫星发射轨道');183.仿真组织函数functionsimu_Satelit直角坐标显示的卫星发射轨道：初始速度越大，规大越大。V不同时，需要的步长和计算步数不一样。19直角坐标显示的卫星发射轨道：初始速度越大，规大越大。193.4线性系统仿真编程前面所介绍的仿真程序是针对(3.1-1)所表示的通用的非线性系统，针对如下所表示的线性状态空间模型(3.4-1)我们虽然也可以用前面介绍的程序仿真，但比较麻烦，实际上系统由矩阵A、B、C、D就可以指定微分方程和输出方程，而不必单独用函数文件来表示。为此，我们单独编写针对系统(3.4-1)的线性系统数值积分法仿真程序。程序仍然由仿真框架和单步积分构成。203.4线性系统仿真编程前面所介绍的仿真程序是针对(3.1-(1)线性系统数值积分仿真框架function[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)%函数功能：用数值积分法仿真对线性系统dx=AX+Bu,y=Cx+Du进行仿真%输入参数：tstart,tstop,h分别是起始时间、结束时间和仿真步长，是标量%x0,u0是状态、输入的初值，都是列向量%A,B,C,D是线性状态空间模型的系数矩阵%InteMethod是数值积分方法，可以是'EUL','RK2','RK4'%ControlFile是控制作用函数，没有控制算法时，可以省略%输出参数：t是仿真结果的时间序列%y是仿真结果系统的输出序列21(1)线性系统数值积分仿真框架function[t,y]=function[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)t=[tstart:h:tstop];%t数一个行序列 cntt=size(t,2);%返回列数cnty=size(C,1);%返回y的维数y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果y0=C*x0+D*u0;%eval([OutputFile,'(tstart,x0,u0)']);%计算初始输出y(:,1)=y0;%将y0作为输出的第1列curx=x0;%当前一步的x curu=u0;%当前一步的u cury=y0;%当前一步的yfori=1:1:cntt-1curu=eval([ControlFile,'(t(i),h,curx,curu,cury,A,B,C,D)']);%计算控制时传递的参数：当前时间，步长，当前状态和输出,模型系数矩阵curx=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%单步积分运算,针对线性模型cury=C*curx+D*curu;%计算输出y(:,i+1)=cury;%将输出加入到输出序列里end22function[t,y]=w_LinearSimu(ts(2)线性系统单步积分函数functionNewX=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%函数功能：单步积分运算,与模型方程有关%输入参数：h是数值积分步长，curx,curu分别是当前的状态和控制向量%InteMethod是积分算法，字符串类型，% 可以取'EUL','RK2','RK4',由于要用于判断，字符串必须等长%A,B是线性状态模型的微分方程的系数矩阵%输出参数：NewX是这一步计算的新的状态向量Bu=B*curu;ifInteMethod=='RK4'k1=A*curx+Bu;k2=A*(curx+0.5*h*k1)+Bu;k3=A*(curx+0.5*h*k2)+Bu;k4=A*(curx+h*k3)+Bu;NewX=curx+h*(k1+2*k2+2*k3+k4)/6;elseifInteMethod=='RK2'k1=A*curx+Bu;%eval([StateModel,'(curt,curx,curu)']);k2=A*(curx+h*k1)+Bu;%eval([StateModel,'(curt+h,curx+h*k1,curu)']);NewX=curx+0.5*h*(k1+k2);else%欧拉法EULNewX=curx+h*(A*curx+Bu);end23(2)线性系统单步积分函数functionNewX=w_L3.5二阶系统传递函数仿真1.模型转换写为状态模型就是基本参数取值为针对本例，线性系统仿真程序进行仿真。243.5二阶系统传递函数仿真1.模型转换写为状态模型就是基本2.仿真组织程序functionsimu_StateSpace%函数功能：对一个状态空间模型进行仿真tstart=0;tstop=20;h=0.1;B=[0;1];C=[1,0];D=[0];x0=[0;0];u0=[1];%单位阶跃输入dampratio=0.1;A=[0,1;-1,-2*dampratio];[t,y1]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,'RK4');dampratio=0.5;A(2,2)=-2*dampratio;[t,y5]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,'RK4');dampratio=1;A(2,2)=-2*dampratio;[t,y10]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,'RK4');dampratio=1.5;A(2,2)=-2*dampratio;[t,y15]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,'RK4');plot(t,y1,'-r',t,y5,':g',t,y10,'*y',t,y15,'.b');%绘制叠加曲线xlabel('x');ylabel('y');title('二阶系统阶跃输入响应');252.仿真组织程序functionsimu_StateSpa二阶系统在阶跃输入作用下，不同阻尼比时的输出曲线。26二阶系统在阶跃输入作用下，不同阻尼比时的输出曲线。263.6面向结构图仿真对下图的系统进行面向结构图的变换，并进行仿真。G1(s)G2(s)G3(s)G4(s)β其中273.6面向结构图仿真对下图的系统进行面向结构图的变换，并进行(1)模型处理每个环节用3个参数来描述，得到的矩阵28(1)模型处理每个环节用3个参数来描述，得到的矩阵28(2)仿真组织程序利用第2章编写的结构图到状态空间模型转换的函数w_bd2ss进行模型处理functionsimu_BlockDiagram%函数功能：面向结构图的仿真G=[2,1,1;0,1,3;5,1,6;0,1,1];W=[0,0,0,0;1,0,-2,0;0,1,0,1;0,0,0,0];W0=[1;0;0;1];[P,Q]=w_bd2ss(G,W,W0);%先通过面向结构图的模型变换得到状态模型tstart=0;tstop=5;h=0.1;C=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];%44个状态都输出D=[0]; x0=[0;0;0;0]; u0=[1];%单位阶跃输入[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,P,Q,C,D,'RK4');%size(y)%plot(t,y(1,:),'-r',t,y(2,:),':g',t,y(3,:),'*m',t,y(4,:),'.b');%绘制叠加曲线plot(t,y(1,:),'-r',t,y(3,:),'*m');%绘制叠加曲线xlabel('time'); ylabel('y'); title('面向结构图的仿真');29(2)仿真组织程序functionsimu_BlockDG=[2,1,1;0,1,3;5,1,6;0,1,1];只绘制y1和y3，可以看出两个输出的区别，它们都趋向于稳定。30G=[2,1,1;0,1,3;5,1,6;0,1,1];只绘同时绘制4个输出时，由于y2,y4变化很大，发散，所以y1和y3看起来好像重叠了。31同时绘制4个输出时，由于y2,y4变化很大，发散，所以y1和附录：plot函数参数取值意义Plot函数的一般调用是：plot(x1,y1,option1,x2,y2,option2,….)32附录：plot函数参数取值意义32第三章Matlab编程(数值积分法仿真)3.1连续系统数值积分法仿真编程思路目的：针对下面的系统编写通用的数值积分法仿真程序(3.1-1)对这样的系统进行仿真，实际上涉及到控制的计算、状态的数值积分计算和输出的计算。3个函数g,f,h确定后，就可以完整地描述一个系统。给定初值：u0，y0,系统中的向量都采用列向量表示33第三章Matlab编程(数值积分法仿真)3.1连续系统数function[t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)%函数功能：用数值积分法仿真%输入参数：tstart,tstop,h分别是起始时间、结束时间和仿真步长，是标量%x0,u0是状态、输入的初值，都是列向量%cnty是输出变量的个数%InteMethod时数值积分方法，可以是'EUL','RK2','RK4'%StateModel是状态模型的文件名%ControlFile是控制作用的文件名%OutputFile是系统输出的文件名%输出参数：t是仿真结果的时间序列%y是仿真结果系统的输出序列(1)数值积分法仿真框架函数34function[t,y]=w_DigiInteSimu(function[t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile,ControlFile)t=[tstart:h:tstop];%t数一个行序列cntt=size(t,2);%返回列数y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果y0=eval([OutputFile,'(tstart,x0,u0)']);%计算初始输出y(:,1)=y0’;%将cury作为输出的第1列curx=x0;%当前一步的xcuru=u0;%当前一步的ucury=y0;%当前一步的yfori=1:1:cntt-1curu=eval([ControlFile,'(t(i),h,curx,curu,cury)']);%计算控制时传递的参数：当前时间，步长，当前状态和输出curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%单步积分运算cury=eval([OutputFile,'(t(i),curx,curu)']);%计算输出y(:,i+1)=cury‘;%将输出加入到输出序列里end35function[t,y]=w_DigiInteSimu(functionNewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel)%函数功能：单步积分运算,与模型方程有关%输入参数：curt是当前时间，h是数值积分步长% curx,curu分别是当前的状态和控制向量%InteMethod是积分算法，字符串类型，可以取'EUL','RK2','RK4'%StateModel是状态模型文件名称，字符串类型%输出参数：NewX是这一步计算的新的状态向量(2)单步数值积分法函数单步数值积分函数只是对微分方程组StateModel进行一步的计算，计算法由InteMethod参数指定，可以上欧拉法，RK2或RK4。36functionNewX=w_StepIntegral(cfunctionNewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel);ifInteMethod=='RK4'k1=eval([StateModel,'(curt,curx,curu)']);k2=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k1,curu)']);k3=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k2,curu)']);k4=eval([StateModel,'(curt+h,curx+h*k3,curu)']);NewX=curx+h*(k1+2*k2+2*k3+k4)/6;elseifInteMethod=='RK2'k1=eval([StateModel,'(curt,curx,curu)']);k2=eval([StateModel,'(curt+h,curx+h*k1,curu)']);NewX=curx+0.5*h*(k1+k2);else%欧拉法EULQk=eval([StateModel,'(curt,curx,curu)']);%eval用来执行一个函数，传递函数名和函数的输入参数NewX=curx+h*Qk;end37functionNewX=w_StepIntegral(c函数w_DigiInteSimu和w_StepIntegral构造了一个数值积分法仿真的框架，并不涉及具体的系统。具体的系统由StateModel,ControlFile,OutputFile参个参数决定，实际上就是三个函数文件名，这三个函数输入输出参数必须遵循特定的格式，在准备好由这3个函数描述的系统后，调用w_DigiInteSimu即可进行仿真。还需要一个调用w_DigiInteSimu进行仿真的程序，它指定模型文件，指定初始参数，并且对仿真结果绘图。38函数w_DigiInteSimu和w_StepIntegra3.2算法稳定性分析仿真编程针对下面的系统，求用解析法、欧拉法和RK4分别求解，计算欧拉法最大允许步长，将步长从0.1逐渐增大，比较三种解的效果。解：1）该系统是稳定的，解析解为2）用欧拉法计算本例时，其步长应该满足3）RK4步长条件式是一个高阶不等式，无法直接求解，只能用试探法确定RK4的步长上限。(3.2-1)393.2算法稳定性分析仿真编程针对下面的系统，求用解析法、欧我们使用3.1介绍的两个基本函数，同时要针对系统模型编写3个函数来描述该系统，最后编写一个实现仿真初值设置、仿真求解和仿真结果显示的函数。functionderX=w_SysStateEqu(curt,curx,curu)%functionderX=w_stateEqu(curt,curx,curu)%函数功能：在此函数里写系统的状态方程%输入参数：curt当前时间，curx和curu是当前状态和控制%输出参数：derX是状态的X的导数值derX(1)=-4*curx(1);（1）状态模型函数w_SysStateEqu在该函数里描述(3.2-1)式的状态微分方程40我们使用3.1介绍的两个基本函数，同时要针对系统模型编写3个functionOutputY=w_SysOutput(curt,curx,curu)%functionOutputY=w_SysOutput(curt,curx,curu)%函数功能：计算系统的输出向量%输入参数：curt是当前时间，curx,curu分别是当前的状态和控制向量%输出参数：OutputY是计算出来的输出向量OutputY(1)=curx(1);%根据具体的模型写输出方程（2）系统输出函数w_SysOutputfunctionControlU=w_SysControl(curt,h,curx,curu,cury)%ControlU=w_sysControl(curt,h,curx,curu,cury)%函数功能：计算系统的控制，如u=PID(t,curx,cury)%输入参数：curt是当前时间，h是数值积分步长％ curx,curu,cury分别是当前的状态、控制和输出向量%输出参数：ControlU是计算出来的控制向量ControlU=curu;%根据实际系统写控制（3）系统控制计算函数w_SysControl41functionOutputY=w_SysOutput(c（4）仿真组织函数functionsimu_Stability%函数功能：稳定性分析仿真tstart=0;tstop=7;h=0.7;StateModel='w_SysStateEqu';ControlFile='w_SysControl';OutputFile='w_SysOutput';x0=[1];u0=[0];cnty=1;[t,y1]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'EUL',StateModel,OutputFile,ControlFile);%用欧拉法求解[t,y4]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,'RK4',StateModel,OutputFile,ControlFile);%用RK4求解y5=exp(-4*t);%计算解析解plot(t,y1,'--m',t,y2,'-.r',t,y5,'-k');%绘制叠加曲线xlabel('时间'); ylabel('y');title('算法稳定性分析仿真');该函数的主要任务是：1)设置仿真起始时间、结束时间、仿真步长、初始的x和u2)设置模型的输出方程、状态方程和控制方程所在函数文件名3)调用w_DigiInteSimu进行仿真计算4)将计算结果绘图进行比较。42（4）仿真组织函数functionsimu_Stabili修改simu_Stability函数里的参数h，就可以作出不同步长时的仿真结果。下面是步长从0.1不断增大时的仿真结果h=0.1时的仿真曲线。1】RK4法的曲线与解析解重合，说明RK4法非常准确。2】欧拉法有较大误差。43修改simu_Stability函数里的参数h，就可以作出不h=0.4时的仿真曲线。1】RK4法还是比较准确。2】欧拉法误差很大，出现衰减震荡。3】两种数值算法的步长条件仍然满足，算法仍然稳定。44h=0.4时的仿真曲线。12h=0.6时的仿真曲线。欧拉法的步长条件已经不满足，仿真解发散h=0.6时的仿真曲线。RK4法的步长条件仍然满足。但误差增大45h=0.6时的仿真曲线。欧拉法的步长条件已经不满足，仿真解发h=0.7时的仿真曲线。RK4法的步长条件已经不满足，仿真结果发散。46h=0.7时的仿真曲线。143.3卫星发射仿真卫星发射运动方程G为重力系数系统是高阶微分方程，不能直接用于仿真计算，需要转化为状态方程。定义状态变量：运动方程用极坐标表示，同时还需要用直角坐标输出。1.模型转换473.3卫星发射仿真卫星发射运动方程G为重力系数系统是高阶微得到系统仿真模型取X-Y平面坐标作为输出状态初值为：v是初始发射速度，可分别取该系统没有控制，主要是研究初始发射速度对轨道的影响。48得到系统仿真模型取X-Y平面坐标作为输出状态初值为：v是初始2.模型描述函数构造functionderX=sat_StateEqu(curt,curx,curu)G=401408;derX(1)=curx(3);derX(2)=curx(4);derX(3)=-G/(curx(1)*curx(1))+curx(1)*curx(4)*curx(4);derX(4)=-2*curx(3)*curx(4)/curx(1);derX=derX’;%转换为列向量(1)状态方程functionOutputY=sat_Output(curt,curx,curu)OutputY(1)=curx(1)*cos(curx(2));OutputY(2)=curx(1)*sin(curx(2));OutputY=OutputY’;%转换为列向量(2)输出方程(3)控制方程由于该系统没有控制，因而采用与3.2节同样的控制函数。492.模型描述函数构造functionderX=sat_S3.仿真组织函数functionsimu_Satelite%函数功能：对卫星发射轨道仿真tstart=0;StateModel='sat_StateEqu';ControlFile='w_SysControl';OutputFile='sat_Output'; u0=[0];cnty=2;h=200;tstop=100*h; v=10;x0=[6400,0,0,v/6400]‘;[t,y10]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解h=200;tstop=50*h; v=9;x0=[6400,0,0,v/6400]’;[t,y9]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解h=100;tstop=60*h; v=8;x0=[6400,0,0,v/6400]‘;[t,y8]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,‘RK4’,StateModel,OutputFile,ControlFile);%用RK4求解plot(y10(1,:),y10(2,:),'-',y9(1,:),y9(2,:),'--',y8(1,:),y8(2,:),':');%绘制叠加曲线xlabel('x'); ylabel('y'); title('卫星发射轨道');503.仿真组织函数functionsimu_Satelit直角坐标显示的卫星发射轨道：初始速度越大，规大越大。V不同时，需要的步长和计算步数不一样。51直角坐标显示的卫星发射轨道：初始速度越大，规大越大。193.4线性系统仿真编程前面所介绍的仿真程序是针对(3.1-1)所表示的通用的非线性系统，针对如下所表示的线性状态空间模型(3.4-1)我们虽然也可以用前面介绍的程序仿真，但比较麻烦，实际上系统由矩阵A、B、C、D就可以指定微分方程和输出方程，而不必单独用函数文件来表示。为此，我们单独编写针对系统(3.4-1)的线性系统数值积分法仿真程序。程序仍然由仿真框架和单步积分构成。523.4线性系统仿真编程前面所介绍的仿真程序是针对(3.1-(1)线性系统数值积分仿真框架function[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)%函数功能：用数值积分法仿真对线性系统dx=AX+Bu,y=Cx+Du进行仿真%输入参数：tstart,tstop,h分别是起始时间、结束时间和仿真步长，是标量%x0,u0是状态、输入的初值，都是列向量%A,B,C,D是线性状态空间模型的系数矩阵%InteMethod是数值积分方法，可以是'EUL','RK2','RK4'%ControlFile是控制作用函数，没有控制算法时，可以省略%输出参数：t是仿真结果的时间序列%y是仿真结果系统的输出序列53(1)线性系统数值积分仿真框架function[t,y]=function[t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,InteMethod,ControlFile)t=[tstart:h:tstop];%t数一个行序列 cntt=size(t,2);%返回列数cnty=size(C,1);%返回y的维数y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果y0=C*x0+D*u0;%eval([OutputFile,'(tstart,x0,u0)']);%计算初始输出y(:,1)=y0;%将y0作为输出的第1列curx=x0;%当前一步的x curu=u0;%当前一步的u cury=y0;%当前一步的yfori=1:1:cntt-1curu=eval([ControlFile,'(t(i),h,curx,curu,cury,A,B,C,D)']);%计算控制时传递的参数：当前时间，步长，当前状态和输出,模型系数矩阵curx=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%单步积分运算,针对线性模型cury=C*curx+D*curu;%计算输出y(:,i+1)=cury;%将输出加入到输出序列里end54function[t,y]=w_LinearSimu(ts(2)线性系统单步积分函数functionNewX=w_LinearStepIntegral(h,curx,curu,InteMethod,A,B);%函数功能：单步积分运算,与模型方程有关%输入参数：h是数值积分步长，curx,curu分别是当前的状态和控制向量%InteMethod是积分算法，字符串类型，% 可以取'EUL','RK2','RK4',由于要用于判断，字符串必须等长%A,B是线性状态模型的微分方程的系数矩阵%输出参数：NewX是这一步计算的新的状态向量Bu=B*curu;ifInteMethod=='RK4'k1=A*curx+Bu;k2=A*(curx+0.5*h*k1)+Bu;k3=A*(curx+0.5*h*k2)+Bu;k4=A*(curx+h*k3)+Bu;NewX=curx+h*(k1+2*k2+2*k3+k4)/6;elseifInteMethod=='RK2'k1=A*curx+Bu;%eval([StateModel,'(curt,curx,curu