第四章连续系统的离散化方法

第四章连续系统的离散化方法第一页，共三十八页，编辑于2023年，星期五其一般公式为称为截断误差例4－1用欧拉法求下述微分方程的数值解。第二页，共三十八页，编辑于2023年，星期五解：因欧拉法的递推公式为所以则递推公式为：若取步长由t=0开始计算可得上式的精确解是数值计算与精确解的比较见表t00.10.20.31.0精确解x(t)10.90909090.83333330.76923070.510.90.8190.57190.4627810第三页，共三十八页，编辑于2023年，星期五4.1.3龙格－库塔法将在点处作台劳级数展开，并取线性部分可得

第四页，共三十八页，编辑于2023年，星期五将

代入式比较各项系数得待定系数个数超过方程个数，必须先设定一个系数，然后即可求得其参数。一般有以下几种取法：1、则一般形式

第五页，共三十八页，编辑于2023年，星期五2、则

一般形式3、则

一般形式以上几种递推公式均称为二阶龙格－库塔公式。是较典型的几个常用算法。其中的方法3又称为预估－校正法，或梯形法。第六页，共三十八页，编辑于2023年，星期五其意义如下：用欧拉法以斜率先求取一点，再由此点求得另一斜率然后，从点开始，既不按该点斜率变化，也不按预估点斜率变化，而是取两者平均值求得校正点，即：。第七页，共三十八页，编辑于2023年，星期五四阶龙格－库塔法的计算公式为：对于用状态方程表示的高阶线性系统

其中状态变量为用四阶龙格－库塔法时，有计算公式为：第八页，共三十八页，编辑于2023年，星期五其中，相应的输出为按上式，取不断递推，

即可求得所需时刻的状态变量和输出值德国学者Felhberg对传统的龙格－库塔法提出了改进，在每一个计算步长内对f()函数进行了六次求值，以保证更高的精度和数值稳定性，假设当前的步长为，定义下面6个变量下一步的状态变量可由下式求出：第九页，共三十八页，编辑于2023年，星期五

四阶/五阶龙格－库塔法系数表016/13525/2161/41/4003/83/329/326656/128251408/256512/131932/2197-7200/21977296/219728561/564302197/41041439/216-83680/513-845/4104-9/50-1/51/2-8/272-3544/25651859/4104-11/402/550这一方法又称为四阶/五阶龙格－库塔法。第十页，共三十八页，编辑于2023年，星期五4.1.4微分方程数值解的MATLAB实现1、该指令适用于一阶常微分方程组如遇到高阶常微分方程，必须先将他们转换成一阶微分方程组，即状态方程方可使用。2、输入参数为定义微分方程组

M－函数文件名，可以在文件名加写@，或用英文格式单引号界定文件名。3、在编辑调试窗口中编写一阶常微分方程组的M－函数文件时，每个微分方程的格式必须与一致，即等号严格以“先自变量t，后函数”的固定顺序输入，表示微分方程的序数。左边为带求函数的一阶导数，右边函数的变量4、输入参数“Tspan”规定了常微分方程的自变量取值范围，它以矩阵[t0,tf]的形式输入，表示自变量5、输入参数x0表示初始条件向量，

第十一页，共三十八页，编辑于2023年，星期五微分方程组中的方程个数必须等于初始条件数，这是求微分方程特解所必须的条件。6、输入参数options表示选项参数（包括tol，trace），可缺省，即取默认值，tol是控制结果精度的选项对ode23()函数取，对ode45()函数取。trace为输出形式控制变量，如果trace不为0，则

会将仿真中间结果逐步地由频幕显示出来，否则将不

显示中间结果7、输出参数[t,x]为微分方程组解函数的列表（t和x都是列矩阵），它包含向量t各节点和与对应向量x的第j个分量值（即第j个方程解），i表示节点序列数。8、输出参数[t,x]缺省时，输出解函数的曲线，即函数及其各的曲线。阶导数求解微分方程的指令还有ode113（多步解法器），ode15s（基于数字微分公式的解法器），ode23s（单步解法器），ode23T（梯形规则的一种自由插值实现），ode23TB（二阶隐式龙格－库塔公式）等。第十二页，共三十八页，编辑于2023年，星期五例4－1求解常微分方程

，初始条件为：解：方法1把二阶微分方程化成两个一阶微分方程组：令

则：

首先编制M文件，并且函数名和M文件名相同。

functionxdot=wffc_1(t,x)%定义输入、输出变量和函数文件名xdot=zeros(2,1);%明确xdot的维数xdot(1)=x(1);%第一个微分方程表示形式xdot(2)=-x(1)+2+t^2/pi;%第二个微分方程表示形式方法2写出系统的状态方程第十三页，共三十八页，编辑于2023年，星期五首先编制M文件，并且函数名和M文件名相同。functionxdot=wffc_1(t,x)%定义t,x,xdot和文件名xdot=[01;-10]*x+[0;1]*(2+t^2/pi);%状态方程的表示形式在命令窗口键入[t,x]=ode45(@wffc_1,[0,10],[-1;1])，可得微分方程的数值解，其前10组数据如下：t=00.01670.03350.05020.06700.15070.23440.31820.40190.5964------x=

-1.00001.0000-0.98281.0501-0.96481.0999-0.94601.1494-0.92631.1986-0.81581.4395-0.68561.6709-0.53631.8917-0.36912.10070.08302.5349------第十四页，共三十八页，编辑于2023年，星期五例4－2求VarderPol微分方程，在初始条件下的数值。解解：将二阶微分方程变换成一阶微分方程组则

第十五页，共三十八页，编辑于2023年，星期五编制M文件，并且函数名和M文件名相同。functionxdot=vdpl(t,x)xdot=zeros(2,1);％赋初值，并规定向量的维数。xdot(1)=x(2);％对第一个微分方程进行描述xdot(2)=2*(1-x(1)^2)*x(2)-x(1);％对第二个微分方程进行描述在命令窗口键入[t,x]=ode45('vdpl',[0,20],[1;1])，可得微分方程的数值解，其前10组数据如下：t=00.05020.10050.15070.20100.31360.42620.53890.65150.7484------x=1.00001.00001.04890.94371.09460.87621.13680.79951.17480.71581.24410.51571.29080.31561.31590.13111.3215-0.02781.3131-0.1431------第十六页，共三十八页，编辑于2023年，星期五若在命令窗口输入ode45(@wffc_1,[0,10],[-1;1])，则可得到系统状态曲线，而不输出数据。如图4－4所示。该系统也可直接用SIMULINK进行仿真。首先建立SIMULINK仿真模型如图4－5所示。设置系统参数如下，将积分器初值设置为系统要求的初值，如图4－6所示。第十七页，共三十八页，编辑于2023年，星期五XY示波器参数设置如图4－7所示。第十八页，共三十八页，编辑于2023年，星期五仿真时间参数设置如图4－8所示。系统仿真曲线如图4－9和4－10所示。第十九页，共三十八页，编辑于2023年，星期五4.2数值解法的稳定性及选择原则4.2.1数值解法的稳定性欧拉法对按欧拉公式计算得这是一个一阶差分方程，Z变换后得，其特征值为则该算法的稳定条件为。算法的稳定域如图4－11所示。

第二十页，共三十八页，编辑于2023年，星期五2.梯形法可见只要原微分方程是稳定的（），应用梯形公式进行数值计算时）。因此梯形公式是恒稳公式。算法的稳定域如图4－12所示必然稳定（。梯形公式也是一种隐式公式，所以隐式算法是一种恒稳算法。第二十一页，共三十八页，编辑于2023年，星期五4.2.2数值解法的选择原则一般来说，选择数值计算方法从以下原则考虑：1、精度问题数值计算的精度主要受下面三类误差影响。截断误差、舍入误差和积累误差。其中截断误差同数值算法的阶次有关，阶次越高，截断误差越小，一般减小步长可减小每一步的截断误差。舍入误差则与计算机字长有关，字长越长，舍入误差越小。积累误差是上述两类误差积累的结果，它同积分时间长短有关。一般积分步长越小，则积累误差越大（在一定的积分时间下）。所以，在一定的数值计算方法下，若从总误差考虑，必定有一个最佳步长值。2、计算速度问题计算速度决定于计算的步数以及每一步积分所需的时间，而每一步的计算时间同数值计算方法有关。例如四阶龙格－库塔法每一步需要求四个斜率值，花费较多的时间（但精度高）。在积分方法确定时，应在保证一定的计算精度下尽量能选用较大的步长（必须满足数值稳定的条件），以缩短积分时间（有时往往选用变步长的算法）。3、数值解的稳定性数值算法实际上就是将微分方程化成差分方程进行求解，对一个稳定的微分方程组，经过变换得到的差分方程不一定稳定。不同的数值计算方法有不同的稳定域。从稳定性角度看，隐式比显式好。第二十二页，共三十八页，编辑于2023年，星期五4.3数值算法中的“病态”问题4.3.1“病态”微分方程例4－3，已知系统的状态方程为其中

，试用MATLAB仿真

解：系统状态方程可写为：可建立系统的SIMULINK仿真模型如图4－13所示。其中第一个积分器输出为，第二个积分器输出为，第三个积分器输出为。第二十三页，共三十八页，编辑于2023年，星期五第二十四页，共三十八页，编辑于2023年，星期五采用四阶龙格－库塔法，选取固定步长的仿真方法，取h=0.01，参数设置由图可见，t=0.1s之后，曲线变化趋于平缓，若仍以h=0.01计算，会使计算时间拖得很长。第二十五页，共三十八页，编辑于2023年，星期五由结果可看出，误差很大，仿真有较大的失真。第二十六页，共三十八页，编辑于2023年，星期五因此，受数值稳定性限制，只能取小步长进行仿真，若系统是一个慢变过程，则计算速度大受影响。究其原因是状态方程的系数矩阵A的特征值差异较大。可求出系统的特征值如下：A=[-2119-20;19-2120;40-40-40]第二十七页，共三十八页，编辑于2023年，星期五A=-2119-2019-212040-40-40>>eig(A)ans=-2.0000-40.0000+40.0000i-40.0000-40.0000i三个特征值为：

而系统的特征值在实际系统中反映了动态过渡过程的作用。

各瞬态分量时间常数第二十八页，共三十八页，编辑于2023年，星期五“病态（stiff）”方程。表述如下：一般线性常微分方程组：的系数矩阵A的特征值具有如下特征

（4－13）则称式（4－13）为病态方程，相应的系统称为病态系统。第二十九页，共三十八页，编辑于2023年，星期五4.3.2“病态”系统的仿真方法采用自动变步长数值计算方法

对于例4－3，

第三十页，共三十八页，编辑于2023年，星期五4.4连续系统状态方程的离散化连续定常系统状态方程为：连续系统状态方程解的一般形式为：若上式中令采样时刻的状态

令

采用零阶保持器，到期间保持器的输出为恒定值。且等于时刻的采样值,并记:

令第三十一页，共三十八页，编辑于2023年，星期五最后得:于是离散化的状态空间方程为：c2d()函数可以立即得出连续系统离散化的模型来，该函数的调用命令为：例4－4已知连续系统的状态方程为：，用MATLAB求T=0.1时，相应的离散状态方程。第三十二页，共三十八页，编辑于2023年，星期五在MATLAB命令窗口输入：>>A=[2-1-1;0-10;021];>>B=[7;2;3];>>C=[124];>>D=0;>>[G,H]=c2d(A,B,0.1)G=1.2214-0.1162-0.116200.9048000.20031.1052H=0.74730.19030.3355

则离散状态方程为：第三十三页，共三十八页，编辑于2023年，星期五MATLAB还提供了c2dm()函数来作类似的变换，其调用格式为：[G,H,Cd,Dd]=c2dm(A,B,C,D,T,‘选项’)它与c2d函数的区别在于，它可以允许用户自己选择变换方法。下面采用两种不同的方法进行离散化：第一种：>>[G,H,Cd,Dd]=c2dm(A,B,C,D,0.1,'zoh')G=1.2214-0.1162-0.116200.9048000.20031.1052H=0.74730.19030.3355Cd=124Dd=0第二种[G,H,Cd,Dd]=c2dm(A,B,C,D,0.1,'tustin')G=1.2222-0.1170-0.117000.9048000.20051.1053H=7.48541.90483.3584Cd=0.11110.22470.4152Dd=1.2364可看出用零阶保持器的离散化得出的结果和直接调用c2d()函数所得的结果相同。第三十四页，共三十八页，编辑于2023年，星期五MATLAB的控制工具箱还给出了离散状态方程转化为连续状态方程的函数d2c()，其调用格式为：

[A,B]=d2c(G,H,T)如对上面采用c2d()得到的结果，采用如下的变换可得到原来的连续状态方程[A,B]=d2c(G,H,0.1)A=2.0000-1.0000-1.00000-1.0000002.00001.0000B=7.00002.00003.0000第三十五页，共三十八页，编辑于2023年，星期五除了d2c()函数外，MATLAB还提供了d2cm()函数，其调用格式为：

[A,B,C,D]=d2cm(G,H,Cd,Dd,T,‘选项’)其中转换方法的选项意义如表4－2所述。对上面采用Tustin变换得到的结果，采用如下的变换可得到原来的连续状态方程[A,B,C,D]=d2cm(G,H,Cd,Dd,0.1,'tustin')A=2.0000-1.0000-1.00000-1.0000002.00001.0000B=7.00002.00003.0000C=1.00002.00004.0000D=-2.2204e-016由转换结果可看出，除了在计算中产生了微小误差外，所得的结果和原连续系统相同。第三十六页，共三十八页，编辑于2023年，星期五下面介绍一种间接的变换方法，首先将要转换的连续传递函数转换成连续的状态方程模型，再对该连续状态方程离散化，然后将相应的离散状态方程再转换成传递函数，即得相应的脉冲传递函数。可编制如下的程序，并命名为：tfc2d.m可实现由连续传递函数表示的系统转换成脉冲传递函数表示的离散系统。function[nd,