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

第四章连续系统的离散化方法4.1常微分方程的数值解法4.1.1数值求解的基本概念已知一个一阶微分方程00(,)()xftxxtx0x10211,,,nntthtthtth12,,,nxxx10tth0()xth应用数值求解的思路是，从初值开始，在一系列时刻，求未知解的近似解，h是计算步长。若x对t的各阶导数都存在。则x(t)在时的解用台劳级数表示为：2(2)()00000()()()()()2!!kkhhxthxthxtxtxtk4.1.2欧拉法000,()ttxtx10tth11()xtx初始条件为计算时的，欧拉法的计算公式为1000(,)xxhftx其一般公式为1(,)kkkkxxhftx0tf0f1fc0t1t23(2)(1)000000(,)(,)(,)2!3!!kknhhhRftxftxftxk称为截断误差例4－1用欧拉法求下述微分方程的数值解。2(0)1xxx解：因欧拉法的递推公式为1(,)kkkkxxhftx所以2(,)kkkftxx则递推公式为：21kkkxxhx若取步长0.1h由t=0开始计算可得2122231010.110.90.90.10.90.8190.8190.10.8190.75190.4627810xxxx上式的精确解是11xt数值计算与精确解的比较见表t00.10.20.31.0精确解x(t)10.90909090.83333330.76923070.510.90.8190.57190.46278104.1.3龙格－库塔法2(2)0000()()()()2!hxthxthxtxt100201021101122(,)(,)()KftxKftbhxbKhxxhaKaK2100000(,)(),2hffxxhftxfttxxtx2K100(,)Kftx将在点处作台劳级数展开，并取线性部分可得20012100(,)ffKftxbhbKhttxxtx将1010020012100(,)[(,)]ffxxahftxahftxbhbhKttxxtx1K2K代入式比较各项系数得122122111,,22aaabab待定系数个数超过方程个数，必须先设定一个系数，然后即可求得其参数。一般有以下几种取法：1、121210,1,2aabb则1021002001(,)(,)22xxhKKftxhhKftxK一般形式12121(,)(,)22kkkkkkxxhKKftxhhKftxK2、1212132,,443aabb则10121002001(3)4(,)22(,)33hxxKKKftxKfthxhK一般形式112121(3)4(,)22(,)33kkkkkkhxxKKKftxKfthxhK3、12121,12aabb则10121002001()2(,)(,)hxxKKKftxKfthxhK一般形式112121()2(,)(,)kkkkkkhxxKKKftxKfthxhK以上几种递推公式均称为二阶龙格－库塔公式。是较典型的几个常用算法。其中的方法3又称为预估－校正法，或梯形法。其意义如下：用欧拉法以斜率先求取一点，再由此点求得另一斜率101xxhK211001(,)(,)KftxfthxhK0x1K2K122KKK10012()2hxxhKxKK然后，从点开始，既不按该点斜率变化，也不按预估点斜率变化，而是取两者平均值求得校正点，即：。0tf0f1f0t1t1001()2hxxff四阶龙格－库塔法的计算公式为：11234(22)6kkhxxKKKK1213243(,)(,)22(,)22(,)kkkkkkkkKftxhhKftxKhhKftxKKfthxhK对于用状态方程表示的高阶线性系统XAXBUYCX0(0)XX其中状态变量为12,,nXxxx用四阶龙格－库塔法时，有(,)ftXAXBU计算公式为：11234(22)6kkhXXKKKK其中，1211322433(,)()(,)()()2222(,)()()2222(,)()()kkkkkkkkkkkkkkkkKftXAXBUthhhhKftXKAXKBUthhhhKftXKAXKBUtKfthXhKAXhKBUth相应的输出为11kkYCX按上式，取0,1,2,,kn即可求得所需时刻11,,,nttt()kXt()kYt的状态变量和输出值hiK德国学者Felhberg对传统的龙格－库塔法提出了改进，在每一个计算步长内对f()函数进行了六次求值，以保证更高的精度和数值稳定性，假设当前的步长为，定义下面6个11(,),1,2,,6;1,2,1iikijjkijKhfxKthiji下一步的状态变量可由下式求出：611kkiiixxK四阶/五阶龙格－库塔法系数表iiji*i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这一方法又称为四阶/五阶龙格－库塔法。4.1.4微分方程数值解的MATLAB实现1、,23('',,0,)txodeFunTspanxoptions,45('',,0,)txodeFunTspanxoptions该指令适用于一阶常微分方程组(,),1,2,iiixftxin如遇到高阶常微分方程，必须先将他们转换成一阶微分方程组，即状态方程方可使用。2、输入参数''Fun为定义微分方程组(,),1,2,iiixftxinM－函数文件名，可以在文件名加写@，或用英文格式单引号界定文件名。3、在编辑调试窗口中编写一阶常微分方程组(,),1,2,iiixftxin的M－函数文件时，每个微分方程的格式必须与(,)iiixftx,itxix1,2,in一致，即等号严格以“先自变量t，后函数”的固定顺序输入，表示微分方程的序数。左边为带求函数的一阶导数，右边函数的变量4、输入参数“Tspan”规定了常微分方程的自变量取值范围，它以矩阵[t0,tf]的形式输入，表示自变量0,tttf5、输入参数x0表示初始条件向量，0(0)(0)(0)xxtxtxt微分方程组中的方程个数必须等于初始条件数，这是求微分方程特解所必须的条件。6、输入参数options表示选项参数（包括tol，trace），可缺省，即取默认值，tol是控制结果精度的选项对ode23()函数取，对ode45()函数取。trace为输出形式控制变量，如果trace不为0，则会将仿真中间结果逐步地由频幕显示出来，否则将不显示中间结果610310itit()jixgt7、输出参数[t,x]为微分方程组解函数的列表（t和x都是列矩阵），它包含向量t各节点和与对应向量x的第j个分量值（即第j个方程解），i表示节点序列数。()xgt()()jjxgt8、输出参数[t,x]缺省时，输出解函数的曲线，即函数及其各的曲线。阶导数求解微分方程的指令还有ode113（多步解法器），ode15s（基于数字微分公式的解法器），ode23s（单步解法器），ode23T（梯形规则的一种自由插值实现），ode23TB（二阶隐式龙格－库塔公式）等。例4－1求解常微分方程220txx，初始条件为：(0)1;(0)1;0,10xxt解：方法1把二阶微分方程化成两个一阶微分方程组：令12,xxxx则：122212xxtxx首先编制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写出系统的状态方程2010(2)101txx首先编制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------22(1)0xxxx(0)1;(0)1xx0,20t例4－2求VarderPol微分方程，在初始条件下的数值。解解：将二阶微分方程变换成一阶微分方程组12,xxxx则12221212*(1)xxxxxx编制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------若在命令窗口输入ode45(@wffc_1,[0,10],[-1;1])，则可得到系统状态曲线，而不输出数据。如图4－4所示。该系统也可直接用SIMULINK进行仿真。首先建立SIMULINK仿真模型如图4－5所示。设置系统参数如下，将积分器初值设置为系统要求的初值，如图4－6所示。XY示波器参数设置如图4－7所示。仿真时间参数设置如图4－8所示。系统仿真曲线如图4－9和4－10所示。4.2数值解法的稳定性及选择原则4.2.1数值解法的稳定性1.欧拉法对0(0)xxxx按欧拉公式计算得1(1),0,1,,1kkkkxxhxhxkn这是一个一阶差分方程，Z变换后得()(1)()zXzhXz，其特征值为1zh11zh11zh2h则该算法的稳定条件为2h。算法的稳定域如图4－11所示。0-1-2Re()hIm()h不稳定域2.梯形法111()()22kkkkkkkhhxxffxxx()(1)(1)()22hhzXzXz1212hzh1212hzh1Re01z可见只要原微分方程是稳定的（），应用梯形公式进行