数值积分法仿真

数值积分法仿真第1页，共63页，2023年，2月20日，星期六Overview数值积分方法的原理是什么？病态系统的特点和仿真算法选取？算法的稳定性分析？第2页，共63页，2023年，2月20日，星期六第一节数字仿真原理在连续系统的仿真中，数值积分法可分为两大类：单步法：以龙格-库塔法为代表多步法：以Adams法为代表数值积分法的要素：基本特性：稳定性空间特性：精度时间特性：速度第3页，共63页，2023年，2月20日，星期六数值积分基本原理连续系统的仿真，主要是对一阶微分方程（组）的求解可见仿真关键是对Qm准确，快速的求解第4页，共63页，2023年，2月20日，星期六步长：将时间t离散t(k)(k=1,2,…n)，相邻两点的距离为步长，即h=t(k+1)-t(k)步进法：数值积分法求近似解根据初始值y0，按照离散的时间序列步进求解。 t0t1t2t3…tn y0y1y2y3…tn计算格式：由y(k)计算出y(k+1)（k=0,1,…,n)的递推公式。数值积分基本名词第5页，共63页，2023年，2月20日，星期六数值积分的基本性能基本性能数值积分算法的性能包含：定性特征：稳定性时间粒度：计算速度空间粒度：计算精度不同的数值积分方法的具有不同的稳定性。同一个模型采用不同的积分算法和不同的积分步长h，稳定性不同。第6页，共63页，2023年，2月20日，星期六计算速度和计算精度各种数值积分方法的差分方程是对原微分方程的近似逼近，并且因为计算机的字长有限，存在明显的截断误差。这些误差都和计算步距h密切相关，所以计算步距是影响计算精度、速度和稳定性的重要因素。h取得较大，计算时间少，截断误差大；h取得较小，截断误差就会减小，但在给定时间范围内，计算次数必然增加，使误差积累增加。第7页，共63页，2023年，2月20日，星期六截断误差、累计舍入误差与步长h截断误差、累计舍入误差与步长h关系如图。图中可知，两种误差对步距的要求是矛盾的，但两者之和有一个最小值，步距最好能选在最小值。然而，实际要做到这一点是很困难的。一般只能根据经验确定一个合理步长区，通常将步长h限制在系统的最小时间常数数量级上。第8页，共63页，2023年，2月20日，星期六引理：泰勒级数:如果f(x)在x0点处任意阶可导,则在该邻域内的n阶泰勒公式为：第二节单步法单步数值积分法的核心就是泰勒级数的近似。第9页，共63页，2023年，2月20日，星期六2.1一阶欧拉法对于一阶微分方程故一般的一阶欧拉法递推形式为：第10页，共63页，2023年，2月20日，星期六一阶欧拉法图示第11页，共63页，2023年，2月20日，星期六2.22阶龙格-库塔对于一阶微分方程第12页，共63页，2023年，2月20日，星期六2阶龙格-库塔第13页，共63页，2023年，2月20日，星期六2阶龙格-库塔第14页，共63页，2023年，2月20日，星期六2阶龙格-库塔故一般的二阶龙格-库塔法递推形式为：2阶龙格-库塔只取到泰勒级数展开式中y的二阶导数项，略去了三阶以上高阶导数项。其截断误差正比于步长h3为纪念提出该方法的德国数学家C.Runge和M.W.Kutta，称这种计算方法为二阶龙格-库塔法。第15页，共63页，2023年，2月20日，星期六2阶龙格-库塔图示第16页，共63页，2023年，2月20日，星期六比较第17页，共63页，2023年，2月20日，星期六高阶龙格-库塔(RK-4)一般在计算精度要求较高的情况下，多使用四阶龙格-库塔法。其计算公式为,其截断误差正比于步长h5第18页，共63页，2023年，2月20日，星期六高阶龙格-库塔(RK-4)第19页，共63页，2023年，2月20日，星期六单步法的特点单步法以上介绍的几种数值积分公式，有一个共同的特点，由于采用了泰勒级数展开，在本次计算中，仅仅用到前一步的计算结果，而不需要利用更前面各步的结果。这类计算方法称为单步法。单步法运算有下列优点：(1)需要存储的数据量少，占用的存储空间少。(2)给定初值，就可启动递推公式进行运算(自启动计算能力)(3)容易实现变步长第20页，共63页，2023年，2月20日，星期六第三节变步长龙格-库塔法步长控制是实现高精度的仿真算法的手段之一。实现步长控制涉及：局部误差估计步长控制策略第21页，共63页，2023年，2月20日，星期六3.1误差估计通常设法寻找一个低一阶的龙格-库塔公式，两者的结果之差可以设为误差。为减少计算量，Ki通常要求公用。Runge-Kutta-Merson法(RK34)第22页，共63页，2023年，2月20日，星期六Runge-Kutta-Fehlberg(RK45)计算公式为5阶6级，误差估计低阶公式为4阶五级，具有四阶误差估计和五阶精度，称为RK45法。RK45被公认为对非病态系统仿真的最有效的方法之一。第23页，共63页，2023年，2月20日，星期六Runge-Kutta-Fehlberg(RK45)iaibijcic*i1016/13525/21621/41/40033/83/329/326656/128251408/2565412/131932/2197-7200/21977296/219728561/564302197/410451439/216-83680/513-847/4104-9/50-1/561/2-8/272-3544/25661859/4104-11/402/550第24页，共63页，2023年，2月20日，星期六RKF-12第25页，共63页，2023年，2月20日，星期六RKS-34（1978，Shamping）第26页，共63页，2023年，2月20日，星期六3.2步长控制步长控制策略一般分为：1）加倍-减半法2）最优步长法第27页，共63页，2023年，2月20日，星期六步长控制：加倍-减半法加倍-减半法第28页，共63页，2023年，2月20日，星期六步长控制：最优步长法第29页，共63页，2023年，2月20日，星期六1.2.2步长控制第30页，共63页，2023年，2月20日，星期六龙格-库塔方法的一般形式各种龙格-库塔法的公式都由两部分组成，一个是上一步结果，另一个是步长乘以各点导数的加权和。平均斜率第31页，共63页，2023年，2月20日，星期六第三节线性多步法单步法运算基于泰勒级数展开法，其特点是：(1)需要存储的数据量少，占用的存储空间少。(2)给定初值(t0,y0)，就可启动递推公式进行运算(自启动计算能力。(3)容易实现变步长积分。可有效平衡计算速度和精度之间的矛盾。多步法的基本原理是多项式拟合利用一个多项式取匹配变量的若干已知值和各阶导数。第32页，共63页，2023年，2月20日，星期六线性多步法原理第33页，共63页，2023年，2月20日，星期六3.1预报公式第34页，共63页，2023年，2月20日，星期六3.1预报公式第35页，共63页，2023年，2月20日，星期六3.1预报公式第36页，共63页，2023年，2月20日，星期六预报举例第37页，共63页，2023年，2月20日，星期六3.2校正公式第38页，共63页，2023年，2月20日，星期六3.2校正公式第39页，共63页，2023年，2月20日，星期六3.2校正公式第40页，共63页，2023年，2月20日，星期六预报-校正举例第41页，共63页，2023年，2月20日，星期六预报-校正举例第42页，共63页，2023年，2月20日，星期六3.3Adams公式根据前面的分析，我们可以将预报和校正公式统一写成：第43页，共63页，2023年，2月20日，星期六显式Adams系数第44页，共63页，2023年，2月20日，星期六隐式Adams系数第45页，共63页，2023年，2月20日，星期六3.4多步法的特点与单步法相比，相同精度下，使用过去多步信息，计算量小。隐式法的精度高，稳定性好，但在计算y(n+k)时需要用到f[y(n+k),t(n+k)]，只能采用迭代法计算。缺点之一是不能自启动，需用单步法计算初始值才能启动计算。第46页，共63页，2023年，2月20日，星期六第四节积分算法的稳定性稳定的系统采用不同的积分算法，其稳定性不同稳定性的测试公式为：当第47页，共63页，2023年，2月20日，星期六4.1一阶Adams法的稳定性分析第48页，共63页，2023年，2月20日，星期六一阶Adams法的稳定性分析第49页，共63页，2023年，2月20日，星期六4.2一般算法的稳定性分析根据上例可得数值积分方法稳定域的一般方法。 设系统测试方程为： 而数值积分公式为： 只有当时，算法稳定。各种数值积分算法的稳定域参见书P96图3.9第50页，共63页，2023年，2月20日，星期六主要算法的稳定性一阶、二阶Admas法为恒稳算法，其他算法条件稳定。除恒稳法外，其他算法的步长h必须限制在最小时间的数量级对龙格-库塔法，阶次k增大，稳定域略微增大。对Admas法，阶次k增大，稳定域反而缩小。第51页，共63页，2023年，2月20日，星期六第五节Matlab实现ODE(OrdinaryDifferentaialequation)解法模型描述算法描述算法仿真第52页，共63页，2023年，2月20日，星期六微分方程模型描述Lorenz曲线filename:mdLorenz.m

functiondx=mdLorenz(t,x)dx=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)];end第53页，共63页，2023年，2月20日，星期六数值积分算法描述matlab中的数值积分算法函数的格式如下:function[tout,yout]=solver(ModelName,tspan,x0,option)第54页，共63页，2023年，2月20日，星期六数值积分算法描述%一阶Euler算法,filename:svEulerfunction[tout,yout]=svEuler(odeFcn,tspan,y0)t0=tspan(1);t1=tspan(2);iflength(tspan)<3,h=(t1-t0)/1000;elseh=tspan(3);tout=[t0:h:t1]';N=length(y0);M=length(tout)-1;tout=[t0:h:t1]';yout=[y0';zeros(M,N)];fori=1:Mk1=h*feval(odeFcn,tout(i),y0);y0=y0+k1;yout(i+1,:)=y0';endend第55页，共63页，2023年，2月20日，星期六数值积分算法描述function[tout,yout]=svRungeKutta4(odeFcn,tspan,y0)t0=tspan(1);t1=tspan(2);iflength(tspan)<3,h=(t1-t0)/1000;elseh=tspan(3);tout=[t0:h:t1]';N=length(y0);M=length(tout)-1;tout=[t0:h:t1]';yout=[y0';zeros(M,N)];fori=1:Mk1=h*feval(odeFcn,tout(i),y0);k2=h*feval(odeFcn,tout(i)+h/2,y0+0.5*k1);k3=h*feval(odeFcn,tout(i)+h/2,y0+0.5*k2);k4=h*feval(odeFcn,tout(i)+h,y0+k3);y0=y0+(k1+2*k2+2*k3+k4)/6;yout(i+1,:)=y0';endend第56页，共63页，2023年，2月20日，星期六微分方程模型描述在Matlab文件中调用方法为:[t,y]=svEuler(@eqLorenz,[0,100],x0);第57页，共63页，2023年，2月20日，星期六ode45ODE45Solvenon-stiffdifferentialequations,mediumordermethod.[TOUT,YOUT]=ODE45(ODEFUN,TSPAN,Y0)withTSPAN=[T0TFINAL]integratesthesystemofdifferentialequationsy'=f(t,y)fromtimeT0toTFINALwithinitialconditionsY0.ODEFUNisafunctionhandle.ForascalarTandavectorY,ODEFUN(T,Y)mustreturnacolumnvectorcorrespondingtof(t,y).EachrowinthesolutionarrayYOUTcorrespondstoatimereturnedinthecolumnvectorTOUT.ToobtainsolutionsatspecifictimesT0,T1,...,TFINAL(allincreasingoralldecreasing),useTSPAN=[T0T1...TFINAL].第58页，共63页，2023年，2月20日，星期六ode45ODE45Solvenon-stiffdifferentialequations,mediumordermethod.

[TOUT,YOUT]=ODE45(ODEFUN,TSPAN,Y0,OPTIONS)solvesasabovewithdefaultintegrationpropertiesreplacedbyvaluesinOPTIONS,anargumentcreatedwiththeODESETfunction.SeeODESETfordetails.Commonlyusedoptionsarescalarrelativeerrortolerance'RelTol'(1e-3bydefault)andvectorofabsoluteerrortolerances'AbsTol'(allcomponents1e-6bydefault).Ifcertaincomponentsofthesolutionmustbenon-negative,useODESETtosetthe'NonNegative'propertytotheindicesofthese第59页，共63页，2023年，2月20日，星期六ode45ODE45Solvenon-stiffdifferentialequations,mediumordermethod.

ODE45cansolveproblemsM(t,y)*y'=f(t,y)withmassmatrixMthatisnonsingular.UseODESETtosetthe'Mass'propertytoafunctionhandleMASSifMASS(T,Y)returnsthevalueofthemassmatrix.Ifthemassmatrixisconstant,thematrixcanbeusedasthevalueofthe'Mass'option.IfthemassmatrixdoesnotdependonthestatevariableYandthefunctionMASSistobecalledwithoneinputargumentT,set'MStateDependence'to'none'.ODE15SandODE23Tcansolveproblemswithsingularmassmatrices.第60页，共63页，2023年，2月20日，星期六ode45

ODE45Solvenon-stiffdifferentialequations,mediumordermethod.Example[t,y]=ode45(@vdp1,[020],[20]);plot(t,y(:,1));solvesthesystemy'=vdp1(t,y),usingthedefaultrelativeerrortolerance1e-3andthedefaultabsolutetoleranceof1e-6foreachcomponent,andplotsthefirstcomponentofthesolution.

ClasssupportforinputsTSPAN,Y0,andtheresultofODEFUN(T,Y):float:doubl