《数值积分法仿真幻灯片.ppt》由会员分享，可在线阅读，更多相关《数值积分法仿真幻灯片.ppt（63页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、数值积分法仿真数值积分法仿真第1页，共63页，编辑于2022年，星期六Overviewu数值积分方法的原理是什么？数值积分方法的原理是什么？u病态系统的特点和仿真算法选取？病态系统的特点和仿真算法选取？u算法的稳定性分析？算法的稳定性分析？第2页，共63页，编辑于2022年，星期六第一节第一节 数字仿真原理数字仿真原理u在连续系统的仿真中，数值积分法可分为两大类：在连续系统的仿真中，数值积分法可分为两大类：单步法：以龙格-库塔法为代表多步法：以Adams法为代表u数值积分法的要素：数值积分法的要素：基本特性：稳定性空间特性：精度时间特性：速度第3页，共63页，编辑于2022年，星期六数值积分基
2、本原理数值积分基本原理u连续系统的仿真，主要是对一阶微分方程（组）的求解连续系统的仿真，主要是对一阶微分方程（组）的求解可见仿真关键是对可见仿真关键是对Qm准确，快速的求解准确，快速的求解第4页，共63页，编辑于2022年，星期六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页，共63页，编辑于2022年
3、，星期六数值积分的基本性能数值积分的基本性能u基本性能基本性能数值积分算法的性能包含：定性特征：稳定性时间粒度：计算速度空间粒度：计算精度不同的数值积分方法的具有不同的稳定性。同一个模型采用不同的积分算法和不同的积分步长h，稳定性不同。第6页，共63页，编辑于2022年，星期六计算速度和计算精度计算速度和计算精度各种数值积分方法的差分方程是对原微分方程的近似逼近，并且因为计算机的字长有限，存在明显的截断误差。这些误差都和计算步距h密切相关，所以计算步距是影响计算精度、速度和稳定性的重要因素。h取得较大，计算时间少，截断误差大；h取得较小，截断误差就会减小，但在给定时间范围内，计算次数必然增加，
4、使误差积累增加。第7页，共63页，编辑于2022年，星期六截断误差、累计舍入误差与步长截断误差、累计舍入误差与步长h截断误差、累计舍入误差与步长h关系如图。图中可知，两种误差对步距的要求是矛盾的，但两者之和有一个最小值，步距最好能选在最小值。然而，实际要做到这一点是很困难的。一般只能根据经验确定一个合理步长区，通常将步长h限制在系统的最小时间常数数量级上。第8页，共63页，编辑于2022年，星期六u引理：引理：泰勒级数:如果f(x)在x0点处任意阶可导,则在该邻域内的n阶泰勒公式为：第二节第二节 单步法单步法u单步数值积分法的核心就是单步数值积分法的核心就是泰勒级数的近似泰勒级数的近似。第9页
5、，共63页，编辑于2022年，星期六2.1 一阶欧拉法一阶欧拉法u对于一阶微分方程对于一阶微分方程u故一般的一阶欧拉法递推形式故一般的一阶欧拉法递推形式为：为：第10页，共63页，编辑于2022年，星期六一阶欧拉法图示一阶欧拉法图示第11页，共63页，编辑于2022年，星期六2.2 2阶龙格阶龙格-库塔库塔u对于一阶微分方程对于一阶微分方程第12页，共63页，编辑于2022年，星期六2阶龙格阶龙格-库塔库塔第13页，共63页，编辑于2022年，星期六2阶龙格阶龙格-库塔库塔第14页，共63页，编辑于2022年，星期六2阶龙格阶龙格-库塔库塔u故一般的二阶龙格故一般的二阶龙格-库塔法递推形式为：
6、库塔法递推形式为：2阶龙格-库塔只取到泰勒级数展开式中y的二阶导数项，略去了三阶以上高阶导数项。其截断误差正比于步长 h3为纪念提出该方法的德国数学家C.Runge和M.W.Kutta，称这种计算方法为二阶龙格-库塔法。第15页，共63页，编辑于2022年，星期六2阶龙格阶龙格-库塔图示库塔图示第16页，共63页，编辑于2022年，星期六比较比较第17页，共63页，编辑于2022年，星期六高阶龙格高阶龙格-库塔库塔(RK-4)u一般在计算精度要求较高的情况下，多使用四阶龙格一般在计算精度要求较高的情况下，多使用四阶龙格-库塔库塔法。其计算公式为法。其计算公式为,其截断其截断误误差正比于步差正比
7、于步长长 h5第18页，共63页，编辑于2022年，星期六高阶龙格高阶龙格-库塔库塔(RK-4)第19页，共63页，编辑于2022年，星期六单步法的特点单步法的特点u单步法单步法以上介绍的几种数值积分公式，有一个共同的特点，由于采用了泰勒级数展开，在本次计算中，仅仅用到前一步的计算结果，而不需要利用更前面各步的结果。这类计算方法称为单步法。u单步法运算有下列优点：单步法运算有下列优点：(1)需要存储的数据量少，占用的存储空间少。(2)给定初值，就可启动递推公式进行运算(自启动计算能力)(3)容易实现变步长第20页，共63页，编辑于2022年，星期六第三节第三节 变步长龙格变步长龙格-库塔法库塔
8、法u步长控制是实现高精度的仿真算法的手段之一。步长控制是实现高精度的仿真算法的手段之一。u实现步长控制涉及：实现步长控制涉及：局部误差估计步长控制策略第21页，共63页，编辑于2022年，星期六3.1 误差估计误差估计u通常设法寻找一个低一阶的龙格通常设法寻找一个低一阶的龙格-库塔公式，两者的结果之库塔公式，两者的结果之差可以设为误差。为减少计算量，差可以设为误差。为减少计算量，Ki通常要求公用。通常要求公用。uRunge-Kutta-Merson法法(RK34)第22页，共63页，编辑于2022年，星期六Runge-Kutta-Fehlberg(RK45)u计算公式为计算公式为5阶阶6级，误
9、差估计低阶公式为级，误差估计低阶公式为4阶五级，具有阶五级，具有四阶误差估计和五阶精度，称为四阶误差估计和五阶精度，称为RK45法。法。uRK45被公认为对非病态系统仿真的最有效的方法之一。被公认为对非病态系统仿真的最有效的方法之一。第23页，共63页，编辑于2022年，星期六Runge-Kutta-Fehlberg(RK45)iaibijcic*i1016/13525/21621/41/40033/83/329/326656/128251408/2565412/131932/2197-7200/21977296/219728561/564302197/410451439/216-83680/
10、513-847/4104-9/50-1/561/2-8/272-3544/25661859/4104-11/402/550第24页，共63页，编辑于2022年，星期六RKF-12第25页，共63页，编辑于2022年，星期六RKS-34（1978，Shamping）第26页，共63页，编辑于2022年，星期六3.2 步长控制步长控制u步长控制策略一般分为：步长控制策略一般分为：1）加倍-减半法2）最优步长法第27页，共63页，编辑于2022年，星期六步长控制：步长控制：加倍-减半法u加倍加倍-减半法减半法第28页，共63页，编辑于2022年，星期六 步长控制：最优步长法步长控制：最优步长法第29
11、页，共63页，编辑于2022年，星期六1.2.2 步长控制步长控制第30页，共63页，编辑于2022年，星期六龙格龙格-库塔方法的一般形式库塔方法的一般形式u各种龙格各种龙格-库塔法的公式都由两部分组成，一个是上一步结果，库塔法的公式都由两部分组成，一个是上一步结果，另一个是步长乘以各点导数的加权和。另一个是步长乘以各点导数的加权和。平均平均斜率斜率第31页，共63页，编辑于2022年，星期六第三节第三节 线性多步法线性多步法u单步法运算基于泰勒级数展开法，其特点是：单步法运算基于泰勒级数展开法，其特点是：(1)需要存储的数据量少，占用的存储空间少。(2)给定初值(t0,y0)，就可启动递推公
12、式进行运算(自启动计算能力。(3)容易实现变步长积分。可有效平衡计算速度和精度之间的矛盾。u多步法的基本原理是多项式拟合多步法的基本原理是多项式拟合利用一个多项式取匹配变量的若干已知值和各阶导数。第32页，共63页，编辑于2022年，星期六线性多步法原理线性多步法原理第33页，共63页，编辑于2022年，星期六3.1 预报公式预报公式第34页，共63页，编辑于2022年，星期六3.1 预报公式预报公式第35页，共63页，编辑于2022年，星期六3.1 预报公式预报公式第36页，共63页，编辑于2022年，星期六预报举例预报举例第37页，共63页，编辑于2022年，星期六3.2 校正公式校正公式
13、第38页，共63页，编辑于2022年，星期六3.2 校正公式校正公式第39页，共63页，编辑于2022年，星期六3.2 校正公式校正公式第40页，共63页，编辑于2022年，星期六预报预报-校正举例校正举例第41页，共63页，编辑于2022年，星期六预报预报-校正举例校正举例第42页，共63页，编辑于2022年，星期六3.3 Adams公式公式u根据前面的分析，我们可以将预报和校正公式统一写成：根据前面的分析，我们可以将预报和校正公式统一写成：第43页，共63页，编辑于2022年，星期六显式显式Adams系数系数第44页，共63页，编辑于2022年，星期六隐式隐式Adams系数系数第45页，共
14、63页，编辑于2022年，星期六3.4 多步法的特点多步法的特点u与单步法相比，相同精度下，使用过去多步信息，计算量与单步法相比，相同精度下，使用过去多步信息，计算量小。小。u隐式法的精度高，稳定性好，但在计算隐式法的精度高，稳定性好，但在计算y(n+k)时需要用到时需要用到fy(n+k),t(n+k)，只能采用迭代法计算。，只能采用迭代法计算。u缺点之一是不能自启动，需用单步法计算初始值才能启动缺点之一是不能自启动，需用单步法计算初始值才能启动计算。计算。第46页，共63页，编辑于2022年，星期六第四节第四节 积分算法的稳定性积分算法的稳定性u稳定的系统采用不同的积分算法，其稳定性不同稳定
15、的系统采用不同的积分算法，其稳定性不同u稳定性的测试公式为：稳定性的测试公式为：u当当第47页，共63页，编辑于2022年，星期六4.1 一阶一阶Adams法的稳定性分析法的稳定性分析第48页，共63页，编辑于2022年，星期六一阶一阶Adams法的稳定性分析法的稳定性分析第49页，共63页，编辑于2022年，星期六4.2 一般算法的稳定性分析一般算法的稳定性分析u根据上例可得数值积分方法稳定域的一般方法。根据上例可得数值积分方法稳定域的一般方法。设系统测试方程为：设系统测试方程为：而数值积分公式为：而数值积分公式为：只有当只有当 时，算法稳定。时，算法稳定。u各种数值积分算法的稳定域参见书各
16、种数值积分算法的稳定域参见书P96图图3.9第50页，共63页，编辑于2022年，星期六主要算法的稳定性主要算法的稳定性u一阶、二阶一阶、二阶Admas法为法为恒稳算法恒稳算法，其他算法，其他算法条件稳定条件稳定。u除恒稳法外，其他算法的步长除恒稳法外，其他算法的步长h必须限制在最小时间的数量级必须限制在最小时间的数量级u对龙格对龙格-库塔法，阶次库塔法，阶次k增大，稳定域略微增大。增大，稳定域略微增大。u对对Admas法，阶次法，阶次k增大，稳定域反而缩小。增大，稳定域反而缩小。第51页，共63页，编辑于2022年，星期六第五节第五节 Matlab实现实现uODE(Ordinary Diff
17、erentaial equation)解法解法模型描述算法描述算法仿真第52页，共63页，编辑于2022年，星期六微分方程模型描述微分方程模型描述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页，共63页，编辑于2022年，星期六数值积分算法描述数值积分算法描述umatlab中的数值积分算法函数的格式如下中的数值积分算法函数的格式如下:function tout,yout=solver(Mo
18、delName,tspan,x0,option)第54页，共63页，编辑于2022年，星期六数值积分算法描述数值积分算法描述%一阶一阶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*fe
19、val(odeFcn,tout(i),y0);y0=y0+k1;yout(i+1,:)=y0;endend 第55页，共63页，编辑于2022年，星期六数值积分算法描述数值积分算法描述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
20、=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页，编辑于2022年，星期六微分方程模型描述微分方程模型描述u在在Matlab文件中调用方法为文件中调用方法为:t,y=svEuler(eqLorenz,0,100,x0);第57页，共63页，编辑于
21、2022年，星期六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,OD
22、EFUN(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页，共63页，编辑于2022年，星期六ode45u ODE45 S
23、olve 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 tole
24、rance 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页，共63页，编辑于2022年，星期六ode45u ODE45 Solve non-stiff differential equ
25、ations,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 opti
26、on.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页，共63页，编辑于2022年，星期六ode45 ODE45 Solve non-stiff differential equations,medium
27、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 OD
28、EFUN(T,Y):float:double,single第61页，共63页，编辑于2022年，星期六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 evalu
29、ating solution:deval ODE examples:rigidode,ballode,orbitode function handles:function_handle 第62页，共63页，编辑于2022年，星期六其他算法简介其他算法简介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页，共63页，编辑于2022年，星期六