[工学]第3章数值积分法在系统仿真中的应用.pptVIP

表3-3 ?四阶/五阶RKF算法系数表 这样，下一步状态变量可以由下式求出: 当然直接采用这一方法是定步长方法，而在Matlab语言中使用的ode45指令采用的是变步长解法， 并引入 （3-32） 误差量来控制步长的大小。 1978年，Shampine提出一套改进的龙格-库塔公式，它每步只计算4次右函数却能够获得4阶精度与3阶误差估计，简称为RKS34算法。具体公式如下： （3-33） 其中： 另外，引入了一个3阶公式 其中： 正好是下一次计算 时的K1，因此只是在第一步要多计算一次右函数f，以后仍每步计算4此右函数f。 RKS34算法的误差估计为： （3-34） ◆ODE解函数： →Ode45：此方法被推荐为首选方法。 →Ode23：这是一个比ode45低阶的方法。 →Ode113：用于更高阶或大的标量计算。 →Ode23t：用于解决难度适中的问题。 →Ode23s：用于解决难度较大的微分方程组。对于系统中存在常量矩阵的情况也有用。 →Ode15s：与ode23相同，但要求的精度更高。 →Ode23tb：用于解决难度较大的问题，对于系统中存在常量矩阵的情况也有用。 值得注意的是，对常微分方程来说，初值问题的数值解法是多种多样的，除了这里介绍的RKF方法外，比较常用的还有Euler法、Adams法、Gear法，它们的侧重应用范围不一样，一些方法侧重于一般问题的仿真，而另一些方法侧重刚性方程的仿真。在Simlink环境中，以内部函数的方式实现了其中一些仿真算法。相关的算法将在后续章节介绍。 1.参数选择函数： ???→ odeset：产生/改变参数结构。 ?? →odeget：得到参数数据。 有许多设置对odeset控制的ODE解是非常有用的，读者可以参见该指令的帮助文件。 2.输出函数： ??→odeplot: 时间列输出函数。 ??→odephas2: 二维相平面输出函数。 ??→odephas3: 三维相平面输出函数。 ??→odeprint: 命令窗打印输出函数。 例：利用ode45求解下面方程组： 这个方程组用在人口动力学中。可以认为是单一化的捕食者---被捕食者模式。例如，狐狸和兔子。x1表示被捕食者，x2表示捕食者。 如果被捕食者有无限的食物，并且不会出现捕食者。于是有 ，这个式子是以指数形式增长的。大量的被捕食者将会使捕食者的数量增长；同样，越来越少的捕食者会使被捕食者 的数量增长。 创建fun函数，将此函数保存在M文件fun.m中： function fun=fun(t,x) fun=[x(1)-0.1*x(1)*x(2)+0.01*t; -x(2)+0.02*x(1)*x(2)+0.04*t]; [t,x]=ode45(‘fun’,[0,20],[30;20]); plot(t,x); xlabel(‘time t0=0,tt=20’); ylabel(‘x values x1(0)=30,x2(0)=20’); grid 然后在Matlab的命令窗口中调用ode45指令或类似指令和画出解的图形： 得到图3-8 图3-8 由函数fun定义的微分方程解的图形 3.2 刚性系统的特点及算法 在工程实践中，在研究化工系统、电子网络、控制系统中，常常会碰见这样的情形 ：一个高阶系统中常由不同的时间常数相互作用着。以惯性导航为例，修正回路时间常数大，稳定回路 时间常数小。导弹、鱼雷等航行器的运动也是如此，质点运动较慢，偏航与俯仰运动较快。所有这些现 象，主要是由于系统模型中的一些小参数，如小时间常数、小质量等存在而引起的。描述这种系统的微 分方程，在数学上常常称为刚性方程，这种系统就称为刚性系统。 现以二阶微分方程组为例，讨论刚性系统以及它在数值求解上的特点。 （3-35） 满足初始条件? 的解为 （3-36） 其中 为A的两个特征值 于是方程(3-36)的解由对应于 和 的分量组成， 但由于 对应系统的最小时间常数，它很快地便无足轻重。表3-4给出系统(3-35)的解方程(3-36)的具体数值。其中 表示忽略最小时间常数 时，系统的简化解。 表3-4 系统(3-35)的解 从表3-4可以看出，系统(3-35)在0.01秒时，系统的精确解和简化解几乎一致，也就是说，在0.01秒以后，我们可以完全忽略 ，精确系统和 简化系统是没有差别的。但在t=0初始时刻，简化系统的初始条件为 ，精确系统和简化系统差别是非常大的，我们又不能够忽略 的影响。这就是系统的边界层效应，即奇异摄动系统。 对于