第四节稳定性,隐格式和其他问题.PDFVIP

第四节：稳定性，隐格式及其他问题 你可能会觉得奇怪，我们为什么要研究如此之多的格式问题。下面的例子会给你一点提示。 ⎧ dx ⎪ dt v 考虑谐振子问题⎨⎪dv ω2 ,不妨取ω=1 。使用简单的2 阶龙格-库塔方法不难作出这个问 − x ⎩dt 题的数值解。我们取简单的两阶和四阶龙格- 库塔格式来做这件事情。取初始条件是 x (0) 0, y (0) 1 ，解就是标准的正弦函数。首先注意到，这个方程的一个特点是总能量 1 2 1 2 x + v 守恒，我们来研究一下数值解的行为。 2 2 当h 很小的时候，所有的数值解看起来都和解析结果一致，例如，取h=0.2 ，二阶龙格-库塔 方法的结果类似于这样： 当我们观察系统能量的时候，发现系统能量在缓慢的改变，这也在预料之中，因为毕竟不是 一个严格解： 然而，当执行足够多的步数后，我们发现，二阶方法得到的解，开始出现类似于振幅指数增 长的现象： 当然，唯一的解释就是直接用我们前面的二阶算法处理这个问题是不合适的。这种算法产生 剧烈振荡或者发散的行为，称为不稳定现象。具体解释这种不稳定行为需要做一点计算，我 们可以用最简单的一阶格式定性理解这个问题，一阶格式是 x x +hy ⎧ n+1 n n ⎨ y n+1 y n −hxn ⎩ ⎧x k r n 我们考虑一个基本的解⎨ n 1 n ，代入后得到r (1±ih) 。 y k r ⎩ n 2 这是什么意思呢？这个解是复数的，我们的计算全是实数。但没有关系，因为任何一个实数 都可以拆成一对共轭复数之和，我们前面写下的基本解，其实是exp(ist) 的一个表达式，而 这个虚指数函数正好就是正弦函数的另外一种表述，换句话说，我们在分析一个振荡模。我 们可以把基本解写成 1+h2 exp[in arctan(h)] 和 1+h2 exp[−in arctan(h)] 。这样，当h 很小的时候，指数上的因子就是±it ，这无疑是正确的，问题是，r 的模是大于1 的。此外， 我们的迭代是线性的，一开始如果某个初始条件被分解成了两个基本解的组合，两个基本解 就会独立发展，总的解永远是它们的组合，因此整个解会变得越来越大。于是，这样的算法 在执行一段时间后必然会发散。 用同样的手法，我们可以解出我们用的二阶和四阶方法的特征根，结果发现，当h 很小的时 候，二阶解法的|r|大于 1，而四阶方法|r|是小于 1 的，也就是说，四阶方法是稳定的，而二 阶方法呈现一个缓慢的发散。 从数值例子可以看到四阶算法产生的能量（因此，振幅）在逐步缩减： 尽管四阶方法是稳定的，但|也带来了某种不太令人满意的副作用。因为这说明随着模拟的 发展，系统的振幅会逐渐减小。这在物理上是很讨厌的，因为它意味着能量不守恒。我们必 须指出的是，通常情况下，这种“耗散”（或者说数值阻尼）并不是有害的，因为它只是消 除某些振荡，只要我们关心的问题不是这些振荡，这就无所谓。不过，现在我们考虑的正好 是一个谐振子问题，物理上它就算振动无数次仍然还是能量守恒的，这耗散就是个很大的麻 烦了。 我们首先考虑一种很古怪的解决问题的方案，就是所谓的“蛙跳”格式。我们简单的把 x 和y 取值的时间点改变一下，也就是说，x 还是取在0,h,2h,…这些时刻点上，而y 则取值在 0.5h,1.5h,….这些“半时刻点”上。稍微考虑一下就知道这时候的简单二阶格式很像前面的 一阶格式，因为二阶精度的推进法是 ( ) ( ) ( / 2) x t h x t v t h ⎧ + + + ⎨ v(t 3/ 2h) v(t h / 2) a(t h) ⎩ + + + + 使用记号xn x (nh), vn v(nh +h / 2) ，于是谐振子问题的迭代格式是 ⎧x n+1 x n +hyn ⎨ y n+1 y n −hxn+1 ⎩ 现在我们用前面的手法分析一下它的稳定性行为，我们发现当h2 的时候，|r|=1，也就是， 在这个格式中，振荡即不增加也