《matlab连续系统模型.ppt》由会员分享，可在线阅读，更多相关《matlab连续系统模型.ppt（60页珍藏版）》请在三一办公上搜索。
1、系统建模与仿真技术第2章 数学模型与系统建模,连续系统模型描述,状态空间模型的概念说明,已知系统如图所示。其状态特性可用下列微分方程描述由上式可知，如果已知uc(t)和i(t)的初始值，以及在t时的外加输入信号，就能够完全唯一地确定在t时的系统状态。将上述微分方程写成矩阵方程的形式，即为状态方程,2.1 连续系统模型描述,通常，用x表示状态矢量，用x1，x2，.表示其分量。对于上式，如令又可写为此处其输出方程为需要指出的是，从理论上讲，描述系统状态的状态变量的选择不是唯一的，可以有无穷多种表示方式。,2.2 MATLAB的连续系统模型,根据上图。在理想条件下，可得到此电路的电压平衡方程式式中，
2、q为电荷量，C为电容。可改写为初始条件为零时，取方程的拉普拉斯变换，取U(s)与Uc(s)之比，即可得到系统的传递函数,2.2 MATLAB的连续系统模型,传递函数模型 在MATLAB中，直接用分子/分母的系数表示num=b0，b1，bm；den=a0，a1，an；sys=tf(num,den)零极点增益模型 z=z0，z1，zm；p=p0，p1，pn；k=k；sys=zpk(z,p,k)状态空间模型 在MATLAB中，该系统可用（A，B，C，D）矩阵组表示。sys=ss(A,B,C,D),2.2 MATLAB的连续系统模型,连续系统数学模型之间的转换ss2tf命令：状态空间模型转换成传递函数
3、模型。num,den=ss2tf(A,B,C,D,iu)式中，iu为输入的序号。ss2zp命令：状态空间模型转换成零极点增益模型。Z,P,K=ss2zp(A,B,C,D,iu)式中，iu为输入的序号。tf2ss命令：传递函数模型转换成状态空间模型。A,B,C,D=tf2ss(num,den)tf2zp命令：将传递函数模型转换成零极点增益模型。Z,P,K=tf2zp(num,den)zp2ss命令：将零极点模型转换成状态空间模型。A,B,C,D=zp2ss(Z,P,K)zp2tf命令：将零极点模型转换成传递函数模型。num,den=zp2tf(Z,P,K),2.2 MATLAB的连续系统模型,机
4、械转动系统如图所示。它由惯性负载和粘性磨擦阻尼器组成。J为转动惯量，f为粘性磨擦系数，为角速度，T为作用到系统上的转矩。对于机械转动系统，其运动方程可写成初始条件为零时，取方程的拉普拉斯变换，取(s)与T(s)之比，即可得到系统的传递函数,2.2 MATLAB的连续系统模型,假设 J=2，f=2.5,k=5num=1;den=2 2.5 5sys=tf(num,den)step(sys);grid onimpulse(sys);grid onpzmap(sys);grid onss(sys)zpk(sys)bode(sys);grid on,ltiview是MATLAB中提供的一个线性系统分析
5、的图形工具。sisotool是一个综合性的用于单输入，单输出系统的分析与设计工具。它为用户设计单输入，单输出系统提供了非常有好的界面。,表示画出sys所描述系统的零极点图,ss用来创建实数或复数的状态空间模型,建立以z为零点,p为极点,k为增益的ZPK模型,2.2 MATLAB的连续系统模型,s+1G(s)=-s3+6s2+11s+6sys=tf(1 1,1 6 11 6)a,b,c,d=tf2ss(1 1,1 6 11 6)sys1=ss(a,b,c,d)T=rot90(eye(size(a)sys11=ss2ss(sys1,T)sys2=minreal(sys)n1=sys2.num1d1
6、=sys2.den1a,b,c,d=tf2ss(n1,d1)sys3=ss(a,b,c,d),MATLAB可以建立多输入-多输出系统,下面是一个两个输入，两个输出的系统模型sys=tf(2 5,2 7 7;3 12 11,1 6 7,1 5 6,1 6 11 6;1 6 11 6,1 5 6)sys(1,1)第1个输入到第1个输出sys(2,1)第1个输入到第2个输出sys(1,2)第2个输入到第1个输出sys(2,2)第2个输入到第2个输出step(sys)bode(sys)step(sys(1,2),2.3 面向结构图的数学模型,在工程实际中，另外一种常用的系统表示方式是系统的结构图（方块
7、图），它是系统中每个元件的功能和信号流号的图解表示。方块图表明了系统中各种元件间的相互关系，能够清楚地表明实际系统中的信号流动情况。在方块图中，通过函数方块，可以将所有的系统变量联系起来。函数方块或简称为方块，是对加到方块上的输入信号的一种运算符号，运算结果以输出量表示。元件的传递函数，通常写进相应的方块中，并以标明信号流向的箭头，将这些方块连接起来。应当指出，信号只能沿箭头方向通过。这样，系统的方块图就清楚地表示了它的单向特性。,2.3 面向结构图的数学模型,如下图表示了一个方块图单元。指向方块的箭头表示输入，而从方块出来的箭头则表示输出。在这些箭头上标明了相应的信号。应当指出，方块输出信号
8、等于输入信号与方块中传递函数的乘积。用方块图表示系统的优点是：只要依据信号的流向，将各元件的方块连结起来，就能够容易地组成整个系统的方块图，通过方块图，还可以评价每一个元件对系统性能的影响。总之，方块图比物理系统本身更容易体现系统的函数功能。方块图包含了与系统动态特性有关的信息，但它不包括与系统物理结构有关的信息。因此，许多完全不同和根本无关的系统，可以用同一个方块图来表示。应当指出，对于一个系统来说，方块图也不是唯一的。由于分析角度的不同，对于同一个系统，可以画出许多不同的方块图。,2.3 面向结构图的数学模型,任何系统，都可以用由方块、相加点和分支点组成的方块图来表示。所谓分支点，就是由方
9、块出来的输出信号，从这一点起同时进入另一个方块或相加点。画方块图的步骤 在绘制系统的方块图时，首先列写描述每一个元件动态特性的方程式。然后假定初始条件等于零，对这些方程式进行拉普拉斯变换，并将每一个拉普拉斯变换方程分别以方块的形式表示出来。最后将这些方块单元结合在一起，以组成完整的方块图。,2.3 面向结构图的数学模型,典型环节的传递函数积分环节惯性环节一阶领先(或迟后)环节比例积分环节二阶振荡环节,2.3 面向结构图的数学模型,2.3 面向结构图的数学模型,2.3 面向结构图的数学模型,对简单系统的建模可直接采用三种基本模型：传递函数、零极点增益、状态空间模型。但实际中经常遇到几个简单系统组
10、合成一个复杂系统。常见形式有并联、串联、闭环及反馈等连接。1 并联 将两个系统按并联方式连接，在MATLAB中可用运算符+实现。命令格式为sys12=sys1+sys2其对应的结果为：Gp(s)=G1(s)+G2(s)2 串联 将两个系统按串联方式连接，在MATLAB中可用运算符*实现。命令格式为sys12=sys1*sys2其对应的结果为：Gs(s)=G1(s)*G2(s),2.3 面向结构图的数学模型,3闭环 将系统通过正负反馈连接成闭环系统，在MATLAB中可用feedback函数实现。命令格式为 numf,denf=feedback(num1,den1,num2,den2,sign)s
11、ign为可选参数，sign=-1为负反馈，而sign=1对应为正反馈。缺省值为负反馈。其对应的结果为4单位反馈 将两个系统按反馈方式连接成闭环系统（对应于单位反馈系统），在MATLAB中可用cloop函数实现。命令格式为numc,denc=cloop(num,den,sign)sign为可选参数，sign=-1为负反馈，而sign=1对应为正反馈。缺省值为负反馈。其对应的结果为,2.3 面向结构图的数学模型,在实际应用中，往往系统是由多个典型模块组成，连接关系也比较复杂。在这种情况下MATLAB的控制系统工具箱中提供了一个脚本文件blkbuild和一个函数connect()来得到具有相互连接关
12、系的结构图的模型。先使用blkbuild得到原始系统的增广状态方程模型（a,b,c,d）,然后输入各个模块的连接关系，建立连接矩阵Q，最后调用connect()函数来得到总的状态方程模型。下面对blkbuild和connect()用法进行说明。,clear,clcn1=1;d1=1 2;n2=1 1;d2=1 2;n3=-1-3;d3=1 4;Q=1 0 3;2 1 0;3 2 0;nblocks=3;blkbuild;%得到原始系统的的增广模型In=1;Out=2A,B,C,D=connect(a,b,c,d,Q,In,Out);sys=tf(ss(A,B,C,D),2.3 面向结构图的数学
13、模型,步骤1，模块排号，首先对结构图中的模块按照信号流动的方向排号步骤2，得到增广矩阵，输入各个模块的信息。具体方法是模块的总个数参数和各模块中的传递函数的参数nblocks=3n1=1;d1=1 0;%注意，变量名不能变。n2=1 1;d2=1 2;n3=-1-3;d3=1 1;blkbuild执行了blkbuild后，在MATLAB的workspace中的得到原始结构图的增广状态方程模型（a,b,c,d）。这实际上是三个输入对三个输出的模型，还不是我们所需要的。步骤3，指定连接关系Q，和输入，输出向量INPUTS，OUTPUTS。Q应有nblocks行，Q的第一列是相应模块的编号，Q的第k
14、行第二列及以后的元素应包含进入k模块的所有信号的信息，与顺序无关。INPUTS由输入信号进入的模块编号构成，OUTPUTS由输出信号流出的模块编号构成。,2.3 面向结构图的数学模型,Q=1 0 3%由第3个模块流出的信号进入第1个模块 2 1 0%由第1个模块流出的信号进入第2个模块 3 2 0;%由第2个模块流出的信号进入第3个模块INPUTS=1%输入端为第1个模块OUTPUTS=2%输出端为第2个模块步骤4，构造整个系统的模型A,B,C,D=connect(a,b,c,d,Q,INPUTS,OUTPUTS)完成整个过程后，得到的是原结构图的状态空间模型。sys=tf(ss(A,B,C,
15、D)%转换为传递函数模型step(sys);grid on,pause%求系统的阶跃响应rlocus(sys);grid on,pause%画系统的根轨迹图bode(sys);grid on%画系统的频率特性图,2.4 面向测量数据的建模方法,对单输入/单输出系统的数学模型，可以用一条曲线来表示。对具有两个输入的系统，其特性可以用一个曲面表示。已知系统的一组数据(通常是通过测量得到)，要得到此系统的模型的函数表达式，称为曲线拟合或曲面拟合。在数据分析上都称为回归分析(Rrgression Analysis)或数据拟合(Data Fitting)。面向测量数据的建模方法常用的有两类，一是插值方法
16、（Interpolation Method）可用于预估在已知数据点中间的函数值，其应用范围相当广泛。二是曲线拟合（Curve Fitting），其目的是利用有限的采样点（Sample Points）来建立数学模型，并通过此模型进行进一步的分析和预测。MATLAB 提供了一系列的内插法和曲线拟合函数，以应付不同的需求。下面的数据是美国1790年至1990年(以10年为单位)的总人口。load census.mat;plot(cdate,pop,o);pause;A=ones(size(cdate),cdate,cdate.2;y=pop;theta=Ayplot(cdate,pop,o,cdat
17、e,A*theta,-),2.4 面向测量数据的建模方法,一维内方法interp1是MATLAB的基本一维内插命令,使用方法是yi=intper1(x,y,xi,method)向量x是数据点的x坐标,向量y是数据点的y坐标,xi是插值点,字符串method则规定插值的方法,共有4种临近点内插法(method=nearest)线性内插法(method=linear)三次样条内插法(method=spline)三次多项式内插法(method=cubic),MATLAB 中interp1使用方法的例子x=0:1:4*pi;y=sin(x).*exp(-x/5);xi=0:0.1:4*pi;y1=in
18、terp1(x,y,xi,nearest);y2=interp1(x,y,xi,linear);y3=interp1(x,y,xi,spline);y4=interp1(x,y,xi,cubic);plot(x,y,o,xi,y1,r,xi,y2,g,xi,y3,b,xi,y4,m);legend(Original,Nearest,Linear,Spline,Cubic);,2.4 面向测量数据的建模方法,已知某平原地区的一条公路经过如下坐标点，请用不同的插值方法绘出这条公路（不考虑公路的宽度）。对于上表给出的数据，请用三种不同插值方法编程绘出这条公路。X(m)0 30 50 70 80 90
19、 120 148 170 180 Y(m)80 64 47 42 48 66 80 120 121 138 X(m)202 212 230 248 268 271 280 290 300 312 Y(m)160 182 200 208 212 210 200 196 188 186 X(m)320 340 360 372 382 390 416 430 478 440 Y(m)200 184 188 200 202 240 246 280 296 308 X(m)420 380 360 340 320 314 280 240 200 Y(m)334 328 334 346 356 360 39
20、2 390 400,2.4 面向测量数据的建模方法,X1=0 30 50 70 80 90 120 148 170 180;Y1=80 64 47 42 48 66 80 120 121 138;X2=202 212 230 248 268 271 280 290 300 312;Y2=160 182 200 208 212 210 200 196 188 186;X3=320 340 360 372 382 390 416 430 478 440;Y3=200 184 188 200 202 240 246 280 296 308;X4=420 380 360 340 320 314 280
21、 240 200;Y4=334 328 334 346 356 360 392 390 400;x=cat(2,X1,X2,X3,X4);y=cat(2,Y1,Y2,Y3,Y4);,t=1:39;dt=1:0.01:39;m,n=size(dt);x1i=interp1(t,x,dt,spline);y1i=interp1(t,y,dt,spline);x2i=interp1(t,x,dt,nearset);y2i=interp1(t,y,dt,nearset);x3i=interp1(t,x,dt,cublic);y3i=interp1(t,y,dt,cublic);subplot(221)
22、plot(x,y,o-);title(折线)subplot(222)plot(x1i,y1i);title(三次样条)subplot(223)plot(x2i,y2i);title(近点插值)subplot(224)plot(x3i,y3i);title(内插),2.4 面向测量数据的建模方法,对于有两个输入的系统的建模问题，可以用二维栅格点内插法MATLAB的interp2用来进行二维栅格点内插,使用方法是zi=intper2(x,y,z,xi,yi,method)其中z是一个矩阵,代表一个函数的高度,矩阵x,y是此函数在方栅格点的x,y坐标,字符串method则规定插值的方法,主要有4种临
23、近点内插法(method=nearest)二维线性内插法(method=bilinear)二维样条内插法(method=spline)二维三次多项式内插法(method=bicubic),2.4 面向测量数据的建模方法,MATLAB的interp2用来进行二维栅格点内插的例子x,y=meshgrid(-3:1:3);z=peaks(x,y);xi,yi=meshgrid(-3:0.25:3);zi1=interp2(x,y,z,xi,yi,nearest);zi2=interp2(x,y,z,xi,yi,bilinear);zi3=interp2(x,y,z,xi,yi,bicubic);su
24、bplot(2,2,1);surf(x,y,z);title(Original);subplot(2,2,2);surf(xi,yi,zi1);title(nearest);subplot(2,2,3);surf(xi,yi,zi2);title(bilinear);subplot(2,2,4);surf(xi,yi,zi3);title(bicubic);,2.4 面向测量数据的建模方法,在某海域测得一些点(x,y)的水深z(单位：米)如下表，水深数据是在 低潮时测得的，船的吃水深度为5米，问在矩形区域（75，200）（50，150）平方米的那些地方要避免进入。x 129.0 140.0 1
25、03.5 88.0 185.5 195.0 105.5 y 7.5 141.5 23.0 147.0 22.5 137.5 85.5 z 4 8 6 8 6 8 8 x 157.5 107.5 77.0 81.0 162.0 162.0 117.5 y-6.5-81.0 3.0 56.5-66.5 84.0-33.5 z 9 9 8 8 9 4 9,2.4 面向测量数据的建模方法,clc;clearx=129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5;y=7.5 141.5 23.0
26、 147.0 22.5 137.5 85.5-6.5-81.0 3.0 56.5-66.5 84.0-33.5;z=4 8 6 8 6 8 8 9 9 8 8 9 4 9;xi,yi=meshgrid(75:5:200,-50:5:200);zi=griddata(x,y,z,xi,yi,v4);mesh(xi,yi,zi);holdmesh(xi,yi,5*ones(size(zi),2.4 面向测量数据的建模方法,三维栅格点内插法，可以对有三个输入变量的系统模型进行分析MATLAB的 interp3 命令可以进行三维栅格点的内插，其一般格式是vi=interp3(x,y,z,v,xi,yi
27、,zi,method)其中，x，y，z，v是三维矩阵，前三者表示数据点的输入部分，v是数据点的输出部分。xi,yi,zi,是内插点，字符串method指定不同的内插方法，共有四种临近点内插（method=nearest）三维线性内插（method=linear）三维样条内插（method=spline）三维三次内插（method=cubic）使用interp3时，矩阵x，y，z必须是严格的递增或递减。一般地说，x，y，z 数据由ndgrid 命令产生，以保证格式的正确性。,2.4 面向测量数据的建模方法,三维栅格点内插法演示x=0:1:9;y=0:1:9;z=0:1:9;%产生0到9之间的10
28、个数x,y,z=meshgrid(x,y,z);%转换成内插的数据格式（网格）v=sin(-x.2-y.2-z.2);slice(x,y,z,v,6,9,2,-2,0);%产生切片图pausexi,yi,zi=meshgrid(0:.5:9,0:.5:9,0:.5:9);%细化切片图vi=interp3(x,y,z,v,xi,yi,zi);slice(xi,yi,zi,vi,6,9,2,-2,0);pause;colorbar%用颜色表示数据的大小,2.4 面向测量数据的建模方法,对于高维系统，MATLAB的 interpn 命令可以进行高维栅格点的内插，其一般格式是vi=interpn(x1
29、,x2,v,y1,y2,method)其中，x1,x2,表示数据点的输入部分，v是数据点的输出部分。y1,y2,是内插点，字符串method指定不同的内插方法，共有四种临近点内插（method=nearest）n 维线性内插（method=linear）n 维样条内插（method=spline）n 维三次内插（method=cubic）使用interpn时，矩阵x，y，z必须是严格的递增或递减。一般地说，x，y，z 数据由ndgrid 命令产生，以保证格式的正确性。,2.4 面向测量数据的建模方法,三角内插法与计算几何MATLAB的提供了一系列的命令解决三角化（Triangulation）及
30、计算几何（Computational Geometry）方面的问题。给定一组平面上的点，delaunay 命令可以返回一组由这些点组成的三角形，且任何一点都不会位于任何一个三角形的外接圆上。load seamount;plot(x,y,.,markersize,12,color,r);hold on;pause;tri=delaunay(x,y);trimesh(tri,x,y,z);hold off;pause;trisurf(tri,x,y,z);axis tight%画出三角内插产生的曲面,2.4 面向测量数据的建模方法,voronoi 图形和delaunay 图形关系非常密切，voro
31、noi 命令可以将数据点所在的平面划分成数个多边形，每个多边形只包含一个数据点。load seamount;voronoi(x,y);MATLAB的convhull 命令可用于计算一组数据点的最小凸多边形（Convex Hulls）load seamount;plot(x,y,.,markersize,12,color,r);hold on;pause;k=convhull(x,y);plot(x(k),y(k);hold off,2.4 面向测量数据的建模方法,线性回归对单输入/单输出系统的数学模型,可以用一条曲线来表示,对具有两个输入的系统,其特性可以用一个曲面表示.已知系统的一组数据(通
32、常是通过测量得到),要得到此系统的模型的函数表达式.可以用回归分析(Rrgression Analysis)的方法。先看一个简单的例子load census.mat;plot(cdate,pop,o);pause;A=ones(size(cdate),cdate,cdate.2;y=pop;theta=Ayplot(cdate,pop,o,cdate,A*theta,-),2.4 面向测量数据的建模方法,下面以peaks函数为例，来说明一般的线性回归的原理。peaks函数的表达式我们假设（1）基底函数已知（2）训练数据包含正态分布的噪声，因此，上述函数可以写成,2.4 面向测量数据的建模方法,
33、其中a,b,c为未知的参数，n为均值为0，方差为1的正态分布噪声。如果要得到30组训练数据，可输入point_n=30;xx,yy,zz=peaks(point_n);subplot(2,2,1);surf(xx,yy,zz);axis tight;pausezz=zz+randn(size(zz);%加入正态分布噪声subplot(2,2,2);surf(xx,yy,zz)axis tight注意两图的差异,2.4 面向测量数据的建模方法,由于噪声很大，所以两图的差异较大，现在我们要用已知的基底函数，找出最佳的参数a,b,cx=xx(:);y=yy(:);z=zz(:);A=(1-x).2.
34、*exp(-x.2-(y+1).2),(x/5-x.3-y.5).*exp(-x.2-y.2),exp(-(x+1).2-y.2);theta=Az%估计参数 pausezz=reshape(A*theta,point_n,point_n);Subplot(2,2,3);surf(xx,yy,zz);axis tight%按照估计的参数绘图,2.4 面向测量数据的建模方法,非线性回归非线性回归分析(Nonlinear Rrgression)是一个比较困难的问题,其原因是1.无法一次得到最优解2.无法保证得到最优解3.必需使用各种非线性的最优化算法4.各种相关的数学性质不明显以一个例子来说明MA
35、TLAB中非用回归的方法。,2.4 面向测量数据的建模方法,假设所用的数学模型是其中a1,a2为线性参数，b1,b2为非线性参数，此模型是非线性的，总平方误差是E是参数a1,a2,b1,b2的非线性函数，我们无法用E对这些参数的导数为零的方法得到E的最小值。推而求其次，我们必须用一般的最优化（Optimization）方法，如：梯度下降法（Gradient Desent）或下山式单纯形搜索法（Downhill Simplex Search）等得到E的最小值。,2.4 面向测量数据的建模方法,假设此函数为fun_e1function E=fun_e1(theta,data)x=data(:,1)
36、;y=data(:,2);model_y=theta(1)*exp(theta(3)*x)+theta(2)*exp(theta(4)*x);E=sum(y-model_y).2);其中theta是参数向量，包含了a1，a2，b1，b2，data是观测得到的数据点，返回的E是总平方误差。要求出E的最小值，可使用MATLAB中的fminsearch命令。计算结果为a1=3.5327 a2=1.2390 b1=-4.3153 b2=-0.6840。,2.4 面向测量数据的建模方法,clc;cleardata=0.1 0.2 0.3 0.4 0.4 0.6 0.7 0.8 0.9 1.0 1.1 1
37、.2 1.3 1.4 1.5 3.5 2.5 1.9 1.8 1.3 1.1 1.0 1.1 0.8 0.6 0.6 0.5 0.4 0.4 0.5;init_theta=0 0 0 0;theta=fminsearch(fun_e1,init_theta,data)%fun_e1的定义见前页x=data(:,1);y=data(:,2);esti_y=theta(1)*exp(theta(3)*x)+theta(2)*exp(theta(4)*x);plot(x,y,ro,x,esti_y,b-);legend(Sample data,Regression curve);,3 经典的连续系统
38、建模方法学,3.1 连续系统离散化原理3.2 龙格库塔法3.3 线性多步法3.4 稳定性分析,3.1 连续系统离散化原理,问题：数字计算机在数值及时间上的离散性，被仿真系统数值及时间上的连续性,连续系统仿真，从本质上说：对原连续系统从时间、数值两个方面进行离散化并选择合适的数值计算方法来近似积分运算离散模型原连续模型？（稳定性、准确性、快速性）,3.1 连续系统离散化原理,对仿真建模方法三个基本要求（1）稳定性：若原连续系统是稳定的，则离散化后得到的仿真模型也应是稳定的。（2）准确性：有不同的评价准则，最基本的准则是 绝对误差准则 相对误差准则（3）快速性：若第n步计算对应的系统时间间隔为hn
39、=tn+1-tn计算机由计算需要的时间为Tn，Tn=hn 称为实时仿真，Tnhn称为超实时仿真 Tnhn 称为亚实时仿真,3.2 龙格-库塔法,微分方程RK2公式RK4公式,3.2 龙格-库塔法,数学模型仿真模型1仿真模型2,3.2 龙格-库塔法,一个高精度的仿真方法必须将步长控制作为手段。实现步长控制涉及到误差估计和步长控制策略两方面的问题。龙格-库塔法的误差估计和步长控制策略的基本思想是：每积分一步都设法估计出本步的积分误差en，然后判断是否满足允许误差E，据此选择相应的步长控制策略。每一步的局部误差通常取以下形式 en=en/(|yn|+1)其中|yn|是利用误差估计式计算出的本步的估计
40、误差。当|yn|较大时，en是相对误差，当|yn|较小时，en 是绝对误差。这样作的目的是避免当y 的值很小时，en变得过大。,3.2 龙格-库塔法,步长控制策略有加倍-减半法：设定一个最小误差限和最大误差限，当估计的局部误差大于最大误差限时，将步长减半，并重新计算这一步，当估计的局部误差介于最小误差和最大误差之间时，步长不变，当估计的局部误差小于最小误差限时，将步长加倍。最优步长法：为了使每个积分步在保证精度的前提下取最大步长（最优步长），可以设法根据本步的误差估计，近似确定下一步的可能最大步长。,3.2 龙格-库塔法,RK方法的误差估计通常是设法找一个低一阶的RK公式，将两个公式计算结果之
41、差作为估计误差。例如Runge-Kutta-Fehlberg法的计算公式是RKF1-2公式用另一个一阶公式来估计误差,3.2 龙格-库塔法,RKM3-4公式：误差估计式用3 阶，计算公式为4 阶。RKM3-4公式误差估计公式,3.2 龙格-库塔法,普通二阶龙格-库塔公式一个计算步内分两子步tn时刻：利用当前的un，yn计算k1,计算一次右端函数需tn+h/2时刻：应计算k2，尽管此时yn+h/2已经得到，但un+1则无法得到。(若对un+1也进行预报加大仿真误差)。仿真执行延迟h/2输出要迟后半个计算步距。,3.2 龙格-库塔法,实时二阶龙格-库塔公式,3.2 龙格-库塔法,3.2 龙格-库塔
42、法,设卫星在空中的运动方程，k=401405，是重力系数。卫星轨道用极坐标表示，研究发射的速度对卫星轨道的影响。将其转化为一阶微分方程组，设 y1=r,y2=则有,3.2 龙格-库塔法,根据卫星的发射速度，可以得到初始条件t=0为y1(0)=6400km,y2(0)=0 y3(0)=0。发射速度分别为 V0=8km/h,10km/h,12km/h。由y4(0)=V0/y1(0)得到 y4(0)=0.00125,0.001653,0.001875.MATLAB的仿真程序function dy=weixing(t,y)dy=y(3);y(4);-401408/y(1)/y(1)+y(1)*y(4)
43、*y(4);-2*y(3)*y(4)/y(1);t,y=ode45(weixing,0 5e4,6400 0 0 0.00125);polar(y(:,2),y(:,1),3.3 线性多步法,龙格库塔法也称为“单步法”，原因是在计算yn+1 时只用到了yn，而不直接使用yn-1，yn-2等项。也就是说，在后一步的计算中，仅仅利用前一步的计算结果。从信息利用的观点来看，单步法对信息的利用率不高。“线性多步法”则是利用多步信息计算下一步的值。线性多步法公式很多，这里介绍显式阿达姆斯法和隐式阿达姆斯法,3.3 线性多步法,线性多步法公式很多，这里介绍显式阿达姆斯法和隐式阿达姆斯法1、阿达姆斯开型公式
44、 yn+1=yn+h55f(xn,yn)59f(xn-1,yn-1)+37f(xn-2,yn-2)9f(xn-3,yn-3)/242、阿达姆斯闭型公式 yn+1=yn+h9f(xn+1,yn+1)+19f(xn,yn)5f(xn-1,yn-1)+f(xn-2,yn-2)/243、阿达姆斯预报-校正公式 y(0)n+1=yn+h55f(xn,yn)59f(xn-1,yn-1)+37f(xn-2,yn-2)9f(xn-3,yn-3)/24yn+1=yn+h9f(xn+1,y(0)n+1)+19f(xn,yn)5f(xn-1,yn-1)+f(xn-2,yn-2)/24,3.4 稳定性分析,微分方程
45、dy/dt=f(t,y)，初始条件y(t0)=y0利用数值积分方法进行仿真时常出现的情况是，本来系统是稳定的，仿真结果却得到不稳定的结论。这种现象通常是因为计算步长过大造成的。因为步长h 选的过大时，计算误差较大，数值积分方法会使各种误差传播出去，引起计算的不稳定。关于计算的稳定性可以这样理解，在实际计算中给出的初值y(t0)=y0不一定很准，同时由于计算机字长有限，在计算中回有舍入误差，特别是有截断误差，所有这些误差都可能在计算中传播出去，引起计算的不稳定。仿真方法选择的基本要求：仿真计算不改变原系统的绝对稳定性。,3.4 稳定性分析,下面以一阶显式阿达姆斯法为例，讨论数值积分方法的稳定域微
46、分方程(测试方程)dy/dt=cy，其中c=a+jb，当 Re(c)0时，原系统是稳定的。用一阶显式阿达姆斯法对它进行仿真yn+1=yn+hf(tn,yn)yn+1=yn+hcyn假设yn（n=0,1,2)为它的仿真解，yn+en（n=0,1,2)为准确解，即(yn+1+en+1)=(yn+en)+hc(yn+en)en+1=en+hcenen+1(1+hc)en=0特征方程为z(1+hc)=0,要求特征方程的根在单位圆内，即|1+hc|1 h2|1/a|也就是说要求积分步长h小于等于系统时间常数的两倍。,3.4 稳定性分析,(测试方程)dy/dt=cy，其中c=a+jb|1+hc|1 我们称它所对应的区域为该算法的稳定域。由上例可以得到确定数值积分方法的稳定域的一般方法。设数值积分公式为 yn+1=p(hc)yn 其中p(hc)是一个关于hc的高次多项式，则当|p(hc)|1时，算法才稳定。根据不同算法的表达式，可以在hc域上画出它们的稳定域。除了隐式一、二阶阿达姆斯法为恒稳定外，其它方法都是条件稳定的。除了恒稳定的方法外，其它方法的步长h 都应限制在系统的最小时间常数的数量级。对龙格库塔法,k大则稳定域略有增加，而对阿达姆斯法来说，k大则稳定域反而缩小。,
