《《仿真教学》ppt课件》由会员分享，可在线阅读，更多相关《《仿真教学》ppt课件（43页珍藏版）》请在金锄头文库上搜索。
1、第三章 连续系统数字仿真 3-1 基于数值积分法的连续系统仿真 3-2 基于离散相似法的连续系统仿真 3-3 系统非线性环节的仿真,3.1 基于数值积分法的连续系统仿真 3.1.1 数值积分基本原理 3.1.2 数值积分方法的选择 3.1.3 基于数值积分法的连续系统仿真 3.1.4 数值积分法的MATLAB函数,连续系统的数学模型一般都能以微分方程的形式给出，所以连续系统仿真算法问题通常可归结为如何用计算机来求解微分方程的问题，也就是对一阶微分方程如何进行求解。数值积分法是解决该问题的重要方法之一。 例如，假设有一系统，它的数学模型可用 （1）式所示微分方程来描述,3.1.1 数值积分基本原
2、理,3.1 基于数值积分的连续系统仿真,三、龙格库塔法（RungeKutta法）,龙格库塔法的基本思想是：用几个点上的函数值的线性组合来代替函数的各阶导数，然后再按泰勒级数展开确定其中的系数。 下面以二阶RK法为例介绍其基本原理。,( 5 ),将K2用二元函数泰勒级数展开式展开，并只取前三 项，则有：,首先假设式（1）的解具有下面的形式：,式中，a1、a2、b1、b2为待定系数。,将 K1、K2 代入式（5），得：,( 6 ),比较式（5）与（6）可以得到：,另一方面，将ym+1在ym附近进行泰勒展开，并只取 前三项，则有：,上述3个方程中有4个未知数，因而解不是唯一的。若限定 ，则可得其中一
3、组解：,将它们代入式（5），可得一组计算公式：,截断误差为 O(h3),截断误差为 O(h5),（6） 梯形公式可看作二阶龙格库塔公式，截断误差正比h3,（2）步长 h可变。,RK的特点：,四、亚当姆斯（Adams）法 单步法在计算ym+1时，只利用前一步的ym的值，经过若干步的计算以后，可以求出一系列的值y1，y2，.，yn。如果充分利用前面多步的值来计算ym+1，则可以达到既提高计算速度又能获得较高精度的目的，这就是多步法的基本思想。在多步法中，应用较广的是Adams法。 Adams法有显式积分和隐式积分两种。 Adams显式一般形式： （5） 式中， 为显式Adams公式系数，部分数据如
4、下表,如二次（Adams）亚当姆斯公式： 三次（Adams）亚当姆斯公式： 即是多步法计算公式。 多步法与单步法相比，欲达相同精度，计算工作量较少，在相同条件下多步法比单步法要快。,计算稳定性,仿真计算时，是否仍然稳定呢？先看下面的例子：,从稳定性理论，我们知道如何去从系统的微分,方程或传递函数去判断该系统的稳定性。那么，对,于一个稳定的连续系统，当用某数值积分方法进行,五、数字仿真中的几个问题,（4）,小结： 一般稳定性与步长关系密切（除恒稳公式之外），若用两种显著不同的步长所得到的数值解有明显差别，则可能是这种数值方法不稳定；反之如果基本相同，则一般视为是稳定的 。,3.1.2 数值积分方
5、法的选择 一、积分方法的选择 1精度问题。 在步长相同的条件下，积分方法的阶数越高，精度越高；另外，多步法的精度比单步法高，隐式算法的精度高于显式算法。因此，当需要高精度时，可采用高阶的多步隐式算法和较小的步长。若精度要求不高，一般可选择低阶算法。 2速度问题。 为加快计算速度，可在精度要求不高时，尽量选择低阶的计算工作量少的方法。,3稳定性。 数值解的稳定性必须保证，否则计算结果将失去真实意义。不同的方法有不同的稳定性，要通过合适选择步长来保证稳定性。,4自启动能力。 单步法有自启动能力，多步法没有自启动能 力，必须借助于单步法启动运算之后，才能开始 工作。一般简单的仿真程序多用单步法。,5
6、步长变化能力。 单步法在整个计算中，步长可在一定范围内变化；而多步法则对步长的变化有严格的要求。若要求仿真时步长可变，最好用单步法。 综上所述，积分方法的选择与多种因素有关，各因素之间又相互影响。究竟选哪一种方法，要由具体系统及具体要求而定。,一般情况下四阶RK法应用最为广泛！,二、积分步长的选择 步长的选择很重要。步长过大会增大截断误差，甚至出现数值不稳定现象， 过小了又因增加了步数，而使舍入误差增大。所以存在一个最佳步长使得总误差最小（如图所示）。 在用经验方法选取步长时，一种方法是根据系统方程中最小时间常数Tmin来决定步长，一般取： h=（0.20.05）Tmin,数值积分计算时，积分
7、步长有固定步长和变步长两种工作方式。,固定步长 就是在整个仿真计算过程中，积分步长 h 始终保持不变。 变步长 就是在仿真积分计算的每一步，根据计算误差的大小改变步长h。 目的：在保持一定的计算精度的前提下，尽可能地选取较大的步长。 方法：首先估计计算误差；判断误差是否在允许的误差0范围内；若在允许误差范围内，则该步计算有效，否则计算无效，改变步长，重新计算。 因此，在变步长的积分计算中，必须解决误差估计方法和步长调整策略两个问题。,三、直接采用MATLAB语言根据算法原理编程,clear t0=0; y0=1/3; h=0.1; N=1.5/h; t(1)=0; y(1)=1/3; for
8、i=2:N+1 k1= -30*y0; y1=y0+h*k1; y0=y1; y(i)=y1; t(i)=t(i-1)+h; end y1,clear t0=0; y0=1/3; h=0.1; N=1.5/h; t(1)=0; y(1)=1/3; for i=2:N+1 k1=-30*y0; k2=-30*(y0+h*k1/2); k3=-30*(y0+h*k2/2); k4=-30*(y0+h*k3); y1=y0+h*(k1+2*k2+2*k3+k4)/6; y0=y1; y(i)=y1; t(i)=t(i-1)+h; end y1,stableEuler.m stableRK.m,RK程
9、序框图,习 题,3-1 用Euler法求初值问题,的数值解。取步长 h=0.1。,Euler法： u（1） =1.7848 RK法： u（1） =1.7321,习 题,3-2 用Euler法（h=0.025）及经典Runge-Kutta法（h=0.1）计算初值问题：,在 t=2 时刻的值，并和精确解比较。,3.1.4 数值积分法的MATLAB函数,MATLAB工具箱提供了各种数值积分方法的常用函数，这些函数的基本功能是用数值计算方法求解常系数微分方程（ODE）或微分方程组。,函数调用格式：,说明： （1）solver为数值积分函数名，即ODE45、ODE23、ODE113、ODE15s、ODE
10、23s。 （2）F是描述常微分方程的函数m文件名。 （3）tspan为一个行向量，描述运算的起止时间。例如 0 20。 （4）y0为初始值。,解决问题的一般步骤： 首先把描述系统的高阶微分方程转变为一阶微分方程组形式。 建立一个函数m文件用于描述该一阶微分方程 组。 再建立一个脚本m文件，调用某个solver函数求解微分方程。,【例3.4-1】 已知二阶微分方程,初始条件：,求时间区间t=0 20微分方程的解。,解：分三个步骤求解 （1）将微分方程表示为一阶微分方程组,（2）建立描述微分方程组的函数m文件vdp.m function dy=vdp(t,y) dy=y(2);(1-y(1)*y(
11、1)*y(2)-y(1);,（3）建立脚本m文件example23.m，调用解题器指令求解y,t,y=ode45(vdp,0 20,0,1); plot(t,y(:,1),r-,t,y(:,2),b:); xlabel(t); ylabel(y); legend(y1,y2);,【例3.4-2】 已知一个刚体在无外力作用下运动方程及初始条件为,初始条件：,求时间区间t=0 12微分方程的解。,解：分两个步骤求解,（1）建立描述微分方程组的函数m文 件rigit.m unction dy=rigit(t,y) dy=y(2)*y(3);-y(1)*y(3);-2*y(1)*y(2);,（2）建立
12、脚本m文件example24.m，调用解题器指令求解y,t,y=ode45(rigit,0 12,0 1 1); plot(t,y(:,1),r-,t,y(:,2),b:,t,y(:,3),g-); xlabel(t); ylabel(y); legend(y1,y2,y3);,几个基本概念：,（1）单步法和多步法 单步法是指计算某时刻数值 ，只需前一时刻 有关信息，它是一种能自启动的算法。多步法是指计算某时刻数值 ，需要 时刻有关信息，它是一种不能自启动的算法。 （2）显式法和隐式法 显式法是指计算 时所需数据均已算出。隐式法是指计算 的算式中含有 时刻的数据。因此在使用隐式公式中，需要用另一公式估计这里未知数据，然后用隐式公式进行迭代，这叫预估校正法。显然这种方法不能字启动。 （3）定步长和变步长 定步长为积分步长 h 在仿真运行过程中始终不变。而积分步长在仿真运行过程中自动修正改变为变步长。,
