数值积分法仿真省公共课一等奖全国赛课获奖课件.pptx

西南交通大学电气工程学院电子信息系西南交通大学电气工程学院电子信息系系统仿真系统仿真System Simulation第 1 页 第 1 页 第三章第三章 数值积分法仿真数值积分法仿真第1页第 2 页 Overviewu数值积分方法原理是什么？数值积分方法原理是什么？u病态系统特点和仿真算法选取？病态系统特点和仿真算法选取？u算法稳定性分析？算法稳定性分析？第2页第 3 页 第一节第一节 数字仿真原理数字仿真原理u在连续系统仿真中，数值积分法可分为两大类：在连续系统仿真中，数值积分法可分为两大类：单步法：以龙格-库塔法为代表多步法：以Adams法为代表u数值积分法要素：数值积分法要素：基本特征：稳定性空间特征：精度时间特征：速度第3页第 4 页 数值积分基本原理数值积分基本原理u连续系统仿真，主要是对一阶微分方程（组）求解连续系统仿真，主要是对一阶微分方程（组）求解可见仿真关键是对可见仿真关键是对Qm准确，快速求解准确，快速求解第4页第 5 页 u步长：步长：将时间t离散t(k)(k=1,2,n)，相邻两点距离为步长，即h=t(k+1)-t(k)u步进法：步进法：数值积分法求近似解依据初始值y0，按照离散时间序列步进求解。t0t1t2t3tny0y1y2y3tnu计算格式：计算格式：由y(k)计算出y(k+1)（k=0,1,n)递推公式。数值积分基本名词数值积分基本名词第5页第 6 页 数值数值积分基本性能积分基本性能u基本性能基本性能数值积分算法性能包含：定性特征：稳定性时间粒度：计算速度空间粒度：计算精度不一样数值积分方法含有不一样稳定性。同一个模型采取不一样积分算法和不一样积分步长h，稳定性不一样。第6页第 7 页 计算速度和计算精度计算速度和计算精度各种数值积分方法差分方程是对原微分方程近似迫近，而且因为计算机字长有限，存在显著截断误差。这些误差都和计算步距h亲密相关，所以计算步距是影响计算精度、速度和稳定性主要原因。h取得较大，计算时间少，截断误差大；h取得较小，截断误差就会减小，但在给定时间范围内，计算次数必定增加，使误差积累增加。第7页第 8 页 截断误差截断误差、累计舍入、累计舍入误差与误差与步长步长h截断误差、累计舍入误差与步长h关系如图。图中可知，两种误差对步距要求是矛盾，但二者之和有一个最小值，步距最好能选在最小值。然而，实际要做到这一点是很困难。普通只能依据经验确定一个合理步长区，通常将步长h限制在系统最小时间常数数量级上。第8页第 9 页 u引理：引理：泰勒级数:假如f(x)在x0点处任意阶可导,则在该邻域内n阶泰勒公式为：第二节第二节 单步法单步法u单步数值积分法关键就是单步数值积分法关键就是泰勒级数近似泰勒级数近似。第9页第 10 页 2.1 一阶欧拉法一阶欧拉法u对于一阶微分方程对于一阶微分方程u故普通一阶欧拉法递推形式故普通一阶欧拉法递推形式为：为：第10页第 11 页 一阶欧拉法图示一阶欧拉法图示第11页第 12 页 2.2 2阶龙格阶龙格-库塔库塔u对于一阶微分方程对于一阶微分方程第12页第 13 页 2阶龙格阶龙格-库塔库塔第13页第 14 页 2阶龙格阶龙格-库塔库塔第14页第 15 页 2阶龙格阶龙格-库塔库塔u故普通二阶龙格故普通二阶龙格-库塔法递推形式为：库塔法递推形式为：2阶龙格-库塔只取到泰勒级数展开式中y二阶导数项，略去了三阶以上高阶导数项。其截断误差正比于步长 h3为纪念提出该方法德国数学家C.Runge和M.W.Kutta，称这种计算方法为二阶龙格-库塔法。第15页第 16 页 2阶龙格阶龙格-库塔图示库塔图示第16页第 17 页 比较比较第17页第 18 页 高阶龙格高阶龙格-库塔库塔(RK-4)u普通在计算精度要求较高情况下，多使用四阶龙格普通在计算精度要求较高情况下，多使用四阶龙格-库塔库塔法。其计算公式为法。其计算公式为,其截断其截断误误差正比于步差正比于步长长 h5第18页第 19 页 高阶龙格高阶龙格-库塔库塔(RK-4)第19页第 20 页 单步法特点单步法特点u单步法单步法以上介绍几个数值积分公式，有一个共同特点，因为采取了泰勒级数展开，在此次计算中，仅仅用到前一步计算结果，而不需要利用更前面各步结果。这类计算方法称为单步法。u单步法运算有以下优点：单步法运算有以下优点：(1)需要存放数据量少，占用存放空间少。(2)给定初值，就可开启递推公式进行运算(自开启计算能力)(3)轻易实现变步长第20页第 21 页 第三节第三节 变步长龙格变步长龙格-库塔法库塔法u步长控制是实现高精度仿真算法伎俩之一。步长控制是实现高精度仿真算法伎俩之一。u实现步长控制包括：实现步长控制包括：局部误差预计步长控制策略第21页第 22 页 3.1 误差预计误差预计u通常设法寻找一个低一阶龙格通常设法寻找一个低一阶龙格-库塔公式，二者结果之差能库塔公式，二者结果之差能够设为误差。为降低计算量，够设为误差。为降低计算量，Ki通常要求公用。通常要求公用。uRunge-Kutta-Merson法法(RK34)第22页第 23 页 Runge-Kutta-Fehlberg(RK45)u计算公式为计算公式为5阶阶6级，误差预计低阶公式为级，误差预计低阶公式为4阶五级，含有阶五级，含有四阶误差预计和五阶精度，称为四阶误差预计和五阶精度，称为RK45法。法。uRK45被公认为对非病态系统仿真最有效方法之一。被公认为对非病态系统仿真最有效方法之一。第23页第 24 页 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页第 25 页 RKF-12第25页第 26 页 RKS-34（1978，Shamping）第26页第 27 页 3.2 步长控制步长控制u步长控制策略普通分为：步长控制策略普通分为：1）加倍-减半法2）最优步长法第27页第 28 页 步长控制：步长控制：加倍-减半法u加倍加倍-减半法减半法第28页第 29 页 步长控制：最优步长法步长控制：最优步长法第29页第 30 页 1.2.2 步长控制步长控制第30页第 31 页 龙格龙格-库塔方法普通形式库塔方法普通形式u各种龙格各种龙格-库塔法公式都由两部分组成，一个是上一步结果，库塔法公式都由两部分组成，一个是上一步结果，另一个是步长乘以各点导数加权和。另一个是步长乘以各点导数加权和。平平均均斜斜率率第31页第 32 页 第三节第三节 线性多步法线性多步法u单步法运算基于泰勒级数展开法，其特点是：单步法运算基于泰勒级数展开法，其特点是：(1)需要存放数据量少，占用存放空间少。(2)给定初值(t0,y0)，就可开启递推公式进行运算(自开启计算能力。(3)轻易实现变步长积分。可有效平衡计算速度和精度之间矛盾。u多步法基本原理是多项式拟合多步法基本原理是多项式拟合利用一个多项式取匹配变量若干已知值和各阶导数。第32页第 33 页 线性多步法原理线性多步法原理第33页第 34 页 3.1 预报公式预报公式第34页第 35 页 3.1 预报公式预报公式第35页第 36 页 3.1 预报公式预报公式第36页第 37 页 预报举例预报举例第37页第 38 页 3.2 校正公式校正公式第38页第 39 页 3.2 校正公式校正公式第39页第 40 页 3.2 校正公式校正公式第40页第 41 页 预报预报-校正举例校正举例第41页第 42 页 预报预报-校正举例校正举例第42页第 43 页 3.3 Adams公式公式u依据前面分析，我们能够将预报和校正公式统一写成：依据前面分析，我们能够将预报和校正公式统一写成：第43页第 44 页 显式显式Adams系数系数第44页第 45 页 隐式隐式Adams系数系数第45页第 46 页 3.4 多步法特点多步法特点u与单步法相比，相同精度下，使用过去多步信息，计算量与单步法相比，相同精度下，使用过去多步信息，计算量小。小。u隐式法精度高，稳定性好，但在计算隐式法精度高，稳定性好，但在计算y(n+k)时需要用到时需要用到fy(n+k),t(n+k)，只能采取迭代法计算。，只能采取迭代法计算。u缺点之一是不能自开启，需用单步法计算初始值才能开启缺点之一是不能自开启，需用单步法计算初始值才能开启计算。计算。第46页第 47 页 第四节第四节 积分算法稳定性积分算法稳定性u稳定系统采取不一样积分算法，其稳定性不一样稳定系统采取不一样积分算法，其稳定性不一样u稳定性测试公式为：稳定性测试公式为：u当当第47页第 48 页 4.1 一阶一阶Adams法稳定性分析法稳定性分析第48页第 49 页 一阶一阶Adams法稳定性分析法稳定性分析第49页第 50 页 4.2 普通算法稳定性分析普通算法稳定性分析u依据上例可得数值积分方法稳定域普通方法。依据上例可得数值积分方法稳定域普通方法。设系统测试方程为：设系统测试方程为：而数值积分公式为：而数值积分公式为：只有当只有当 时，算法稳定。时，算法稳定。u各种数值积分算法稳定域参见书各种数值积分算法稳定域参见书P96图图3.9第50页第 51 页 主要算法稳定性主要算法稳定性u一阶、二阶一阶、二阶Admas法为法为恒稳算法恒稳算法，其它算法，其它算法条件稳定条件稳定。u除恒稳法外，其它算法步长除恒稳法外，其它算法步长h必须限制在最小时间数量级必须限制在最小时间数量级u对龙格对龙格-库塔法，阶次库塔法，阶次k增大，稳定域略微增大。增大，稳定域略微增大。u对对Admas法，阶次法，阶次k增大，稳定域反而缩小。增大，稳定域反而缩小。第51页第 52 页 第五节第五节 Matlab实现实现uODE(Ordinary Differentaial equation)解法解法模型描述算法描述算法仿真第52页第 53 页 微分方程模型描述微分方程模型描述uLorenz曲线曲线filename:mdLorenz.m function dx=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页第 54 页 数值积分算法描述数值积分算法描述umatlab中中数值积分算法函数数值积分算法函数格式以下格式以下:function tout,yout=solver(ModelName,tspan,x0,option)第54页第 55 页 数值积分算法描述数值积分算法描述%一阶一阶Euler算法算法,filename:svEulerfunction tout,yout=svEuler(odeFcn,tspan,y0)t0=tspan(1);t1=tspan(2);if length(tspan)3,h=(t1-t0)/1000;else h=tspan(3);tout=t0:h:t1;N=length(y0);M=length(tout)-1;tout=t0:h:t1;yout=y0;zeros(M,N);for i=1:Mk1=h*feval(odeFcn,tout(i),y0);y0=y0+k1;yout(i+1,:)=y0;endend 第55页第 56 页 数值积分算法描述数值积分算法描述function tout,yout=svRungeKutta4(odeFcn,tspan,y0)t0=tspan(1);t1=tspan(2);if length(tspan)3,h=(t1-t0)/1000;else h=tspan(3);tout=t0:h:t1;N=length(y0);M=length(tout)-1;tout=t0:h:t1;yout=y0;zeros(M,N);for i=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页第 57 页 微分方程模型描述微分方程模型描述u在在Matlab文件中调用方法为文件中调用方法为:t,y=svEuler(eqLorenz,0,100,x0);第57页第 58 页 ode45u ODE45 Solve non-stiff differential equations,medium order method.TOUT,YOUT=ODE45(ODEFUN,TSPAN,Y0)with TSPAN=T0 TFINAL integrates the system of differential equations y=f(t,y)from time T0 to TFINAL with initial conditions Y0.ODEFUN is a function handle.For a scalar T and a vector Y,ODEFUN(T,Y)must return a column vector corresponding to f(t,y).Each row in the solution array YOUT corresponds to a time returned in the column vector TOUT.To obtain solutions at specific times T0,T1,.,TFINAL(all increasing or all decreasing),use TSPAN=T0 T1.TFINAL.第58页第 59 页 ode45u ODE45 Solve non-stiff differential equations,medium order method.TOUT,YOUT=ODE45(ODEFUN,TSPAN,Y0,OPTIONS)solves as above with default integration properties replaced by values in OPTIONS,an argument created with the ODESET function.See ODESET for details.Commonly used options are scalar relative error tolerance RelTol(1e-3 by default)and vector of absolute error tolerances AbsTol(all components 1e-6 by default).If certain components of the solution must be non-negative,use ODESET to set the NonNegative property to the indices of these第59页第 60 页 ode45u ODE45 Solve non-stiff differential equations,medium order method.ODE45 can solve problems M(t,y)*y=f(t,y)with mass matrix M that is nonsingular.Use ODESET to set the Mass property to a function handle MASS if MASS(T,Y)returns the value of the mass matrix.If the mass matrix is constant,the matrix can be used as the value of the Mass option.If the mass matrix does not depend on the state variable Y and the function MASS is to be called with one input argument T,set MStateDependence to none.ODE15S and ODE23T can solve problems with singular mass matrices.第60页第 61 页 ode45 ODE45 Solve non-stiff differential equations,medium order method.Example t,y=ode45(vdp1,0 20,2 0);plot(t,y(:,1);solves the system y=vdp1(t,y),using the default relative error tolerance 1e-3 and the default absolute tolerance of 1e-6 for each component,and plots the first component of the solution.Class support for inputs TSPAN,Y0,and the result of ODEFUN(T,Y):float:double,single第61页第 62 页 ode45 ODE45 Solve non-stiff differential equations,medium order method.See also other ODE solvers:ode23,ode113,ode15s,ode23s,ode23t,ode23tb implicit ODEs:ode15i options handling:odeset,odeget output functions:odeplot,odephas2,odephas3,odeprint evaluating solution:deval ODE examples:rigidode,ballode,orbitode function handles:function_handle 第62页第 63 页 其它算法介绍其它算法介绍u ODE23 Solve non-stiff differential equations,low order method.uODE15I Solve fully implicit differential equations,variable order method.uODE23T Solve moderately stiff ODEs and DAEs,trapezoidal rule.u ODE23S Solve stiff differential equations,low order method.uODE23TB Solve stiff differential equations,low order method.第63页