第六章 常微分方程2014.ppt

* * 第8 节 数值稳定性 由上节的讨论可知，步长大小的选择会直接影响到局部截断误差，本节则要介绍步长大小的选择对数值计算方法稳定性的影响。 数值计算方法稳定性：任一步产生的误差在以后均能逐步衰减，则称这种方法是稳定的。 对于稳定的数值计算方法，积分步长的选择将只决定于局部截断误差。 为了分析步长大小的选择对数值计算方法稳定性的影响，引入下面简单的测试方程： 容易得出其解为 当λ0时，x(t)随着时间的推移而趋向于零；当λ0时，x(t)随着时间的推移而趋向无穷大。下面将讨论应用欧拉法、后退欧拉法、多步法数值积分计算时的稳定区域。 对测试方程应用欧拉法： 有 可以看出只有当|1+hλ|1时，才能使当λ0时，x(t)随着时间的推移而趋向于零。 所以对于λ0，系统稳定的条件是hλ落在以（－1，0）为圆心的单位圆内，如图所示。 可以得出λ值越大，步长就必须选择得越小。 相似地，对测试方程使用后退欧拉法： 有 因为当λ0时，x(t)随着时间的推移而趋向于零，所以必须使|1－hλ| 1。即对于λ0，系统稳定的条件是hλ落在以（1，0）为圆心的单位圆外，如图所示。 可以看出使用后退欧拉法时积分步长可以取得任意大，而不至于影响到解的稳定性，这时步长的选择只取决于局部截断误差。 对测试方程应用多步法： 整理上式得 相应的特征方程为 设z1,z2,…,zp+1为上式的根，则有： 当λ0时，xn+1应随着n趋向于无穷而趋向于零，只有当|zj|1（j=1,2,…,p+1）时才能满足。 定义绝对稳定性为整体误差随着n的增大而减少。称多步法是绝对稳定的，如果对于一给定的hλ，方程的根满足|zi|1（i=1,2,…,k）。下面将推导多步法的稳定区域。 令 其中 而 由于z为复数，其可表示为 最终有 多步法的绝对稳定区域的边界可由上式的函数在复平面中绘出，其中hλ(θ)是从0到2π。 绘出三阶的Gear法和三阶Adam法（隐式和显式）的绝对稳定区域。 解： Gear法：令p=k-1，b0=b1=…=0代入上式得： 将三阶的Gear法的系数代入得： 将θ从0取到2π，三阶的Gear法的绝对稳定域即可绘出，如图中阴影区域。 三阶的Gear法的绝对稳定域 令p=k-1，a1=a2=…=ap，得： 同样，代入三阶系数有： 所以三阶的Adam’s-moulton法的绝对稳定域如图中阴影区域所示。 三阶的Adam’s-moulton法的绝对稳定域 令p=k-1，b-1＝0，a1=a2=…=ap=0，同理有 三阶的Adam’s-Bashforth法的绝对稳定域如图中阴影区域所示。 三阶的Adam’s-Bashforth法的绝对稳定域 本例用图像表明了隐式法和显式法的主要差别，隐式法（Gear法和Adam’s-moulton法）的绝对稳定区域要比显式法（Adam’s-Bashforth法）大得多。如Gear法的绝对稳定区域几乎包括整个hλ的左半平面，所以积分步长可以取得任意大，而不至于影响到解的稳定性。正因为这个原因，各种商业软件中往往会采用隐式法，其积分步长的选择只取决于局部截断误差。 另外，随着阶数的增加，绝对稳定区域会不断缩小，但精确度却会增加。 第9 节 刚性系统 刚性系统是指这样一个系统，其时域中的动态特性有变化得非常快的部分，也有变化得很慢的部分。用数学模型表示，即系统的特征值相差很大。 如下式表示的系统即为刚性系统： 特征矩阵为 则相应的特征值为λ1=-1,λ2=-50。 其解为 如图所示，每个状态变量的响应都有快变和慢变部分，快变部分主要是初始的响应，而慢变部分则表现为之后的动态特性。 快变和慢变部分 当考虑用数值方法求解时，由于特征值相差太大，会给步长的选择造成困难。如使用欧拉法求解时，为使求解过程绝对稳定，则要使|1+hλ1|1和|1+hλ2|1， 即要满足0 h2，0 h0.04。所以h要选得小于0.04，这会导致在计算慢变部分时步数过多，效率低下（当λ更大时，可能会导致计算时间长得无法接受)。 所以为了有效而精确地求解刚性微分方程，需要选择合适的算法。 试对下面刚性系统分别使用三阶Adam’s-Bashforth法、Adam’s-moulton法、Gear法，并对结果进行比较。 解：上面曾求出该系统的解为 计算实例 选用h=0.0111s，计算结果如图所示。 可以看出，尽管步