数值积分法仿真公开课一等奖优质课大赛微课获奖课件

1、第三章 数值积分法仿真第1页Overview数值积分方法原理是什么？病态系统特点和仿真算法选取？算法稳定性分析？第2页第一节 数字仿真原理在连续系统仿真中，数值积分法可分为两大类：单步法：以龙格-库塔法为代表多步法：以Adams法为代表数值积分法要素：基本特征：稳定性空间特征：精度时间特征：速度第3页数值积分基本原理连续系统仿真，主要是对一阶微分方程（组）求解可见仿真关键是对Qm准确，快速求解第4页步长：将时间t离散t(k)(k=1,2,n)，相邻两点距离为步长，即h=t(k+1)-t(k)步进法：数值积分法求近似解依据初始值y0，按照离散时间序列步进求解。t0t1t2t3tny0y1y2y3

2、tn计算格式：由y(k)计算出y(k+1)（k=0,1,n)递推公式。数值积分基本名词第5页数值积分基本性能基本性能数值积分算法性能包含：定性特征：稳定性时间粒度：计算速度空间粒度：计算精度不一样数值积分方法含有不一样稳定性。同一个模型采取不一样积分算法和不一样积分步长h，稳定性不一样。第6页计算速度和计算精度各种数值积分方法差分方程是对原微分方程近似迫近，而且因为计算机字长有限，存在显著截断误差。这些误差都和计算步距h亲密相关，所以计算步距是影响计算精度、速度和稳定性主要原因。h取得较大，计算时间少，截断误差大；h取得较小，截断误差就会减小，但在给定时间范围内，计算次数必定增加，使误差积累增

3、加。第7页截断误差、累计舍入误差与步长h截断误差、累计舍入误差与步长h关系如图。图中可知，两种误差对步距要求是矛盾，但二者之和有一个最小值，步距最好能选在最小值。然而，实际要做到这一点是很困难。普通只能依据经验确定一个合理步长区，通常将步长h限制在系统最小时间常数数量级上。第8页引理：泰勒级数:假如f(x)在x0点处任意阶可导,则在该邻域内n阶泰勒公式为：第二节 单步法单步数值积分法关键就是泰勒级数近似。第9页2.1 一阶欧拉法对于一阶微分方程故普通一阶欧拉法递推形式为：第10页一阶欧拉法图示第11页2.2 2阶龙格-库塔对于一阶微分方程第12页2阶龙格-库塔第13页2阶龙格-库塔第14页2阶

4、龙格-库塔故普通二阶龙格-库塔法递推形式为：2阶龙格-库塔只取到泰勒级数展开式中y二阶导数项，略去了三阶以上高阶导数项。其截断误差正比于步长 h3为纪念提出该方法德国数学家C. Runge和M. W .Kutta，称这种计算方法为二阶龙格-库塔法。第15页2阶龙格-库塔图示第16页比较第17页高阶龙格-库塔(RK-4)普通在计算精度要求较高情况下，多使用四阶龙格-库塔法。其计算公式为,其截断误差正比于步长 h5第18页高阶龙格-库塔(RK-4)第19页单步法特点单步法以上介绍几个数值积分公式，有一个共同特点，因为采取了泰勒级数展开，在此次计算中，仅仅用到前一步计算结果，而不需要利用更前面各步结

5、果。这类计算方法称为单步法。单步法运算有以下优点：(1) 需要存放数据量少，占用存放空间少。(2)给定初值，就可开启递推公式进行运算(自开启计算能力)(3) 轻易实现变步长第20页第三节 变步长龙格-库塔法步长控制是实现高精度仿真算法伎俩之一。实现步长控制包括：局部误差预计步长控制策略第21页3.1 误差预计通常设法寻找一个低一阶龙格-库塔公式，二者结果之差能够设为误差。为降低计算量，Ki通常要求公用。Runge-Kutta-Merson法 (RK34)第22页Runge-Kutta-Fehlberg(RK45)计算公式为5阶6级，误差预计低阶公式为4阶五级，含有四阶误差预计和五阶精度，称为R

6、K45法。RK45被公认为对非病态系统仿真最有效方法之一。第23页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页RKF-12第25页RKS-34（1978，Shamping）第26页3.2 步

7、长控制步长控制策略普通分为：1）加倍-减半法2）最优步长法第27页步长控制：加倍-减半法加倍-减半法第28页 步长控制：最优步长法第29页1.2.2 步长控制第30页龙格-库塔方法普通形式各种龙格-库塔法公式都由两部分组成，一个是上一步结果，另一个是步长乘以各点导数加权和。平均斜率第31页第三节 线性多步法单步法运算基于泰勒级数展开法，其特点是：(1) 需要存放数据量少，占用存放空间少。(2)给定初值(t0,y0)，就可开启递推公式进行运算(自开启计算能力。(3) 轻易实现变步长积分。可有效平衡计算速度和精度之间矛盾。多步法基本原理是多项式拟合利用一个多项式取匹配变量若干已知值和各阶导数。第3

8、2页线性多步法原理第33页3.1 预报公式第34页3.1 预报公式第35页3.1 预报公式第36页预报举例第37页3.2 校正公式第38页3.2 校正公式第39页3.2 校正公式第40页预报-校正举例第41页预报-校正举例第42页3.3 Adams公式依据前面分析，我们能够将预报和校正公式统一写成：第43页显式Adams系数第44页隐式Adams系数第45页3.4 多步法特点与单步法相比，相同精度下，使用过去多步信息，计算量小。隐式法精度高，稳定性好，但在计算y(n+k)时需要用到fy(n+k),t(n+k)，只能采取迭代法计算。缺点之一是不能自开启，需用单步法计算初始值才能开启计算。第46页

9、第四节 积分算法稳定性稳定系统采取不一样积分算法，其稳定性不一样稳定性测试公式为：当第47页4.1 一阶Adams法稳定性分析第48页一阶Adams法稳定性分析第49页4.2 普通算法稳定性分析依据上例可得数值积分方法稳定域普通方法。设系统测试方程为：而数值积分公式为：只有当 时，算法稳定。各种数值积分算法稳定域参见书P96图3.9第50页主要算法稳定性一阶、二阶Admas法为恒稳算法，其它算法条件稳定。除恒稳法外，其它算法步长h必须限制在最小时间数量级对龙格-库塔法，阶次k增大，稳定域略微增大。对Admas法，阶次k增大，稳定域反而缩小。第51页第五节 Matlab实现ODE (Ordina

10、ry Differentaial equation)解法模型描述算法描述算法仿真第52页微分方程模型描述Lorenz曲线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页数值积分算法描述matlab中数值积分算法函数格式以下:function tout,yout =solver(ModelName,tspan,x0,option)第54页数值积分算法描述%一阶Euler算法,filename:svEule

11、rfunction 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页数值积分算法描述function tou

12、t,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

13、,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页微分方程模型描述在Matlab文件中调用方法为:t,y=svEuler(eqLorenz,0,100,x0);第57页ode45 ODE45 Solve non-stiff differential equations, medium order method. TOUT,YOUT = ODE45(ODEFUN,TSPAN,Y0) with TSPAN = T0 TFI

14、NAL 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 ret

15、urned 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页ode45 ODE45 Solve non-stiff differential equations, medium order method. TOUT,YOUT = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default i

16、ntegration 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 co

17、mponents of the solution must be non-negative, use ODESET to set the NonNegative property to the indices of these第59页ode45 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

18、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, se

19、t MStateDependence to none. ODE15S and ODE23T can solve problems with singular mass matrices. 第60页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页ode45 ODE45 Solve non-stiff differenti