16.901 讲义笔记 2002.2.11 主题： * 简化为一阶系统 * 线性方程组系统的特性 * 刚度 * 隐含模型与显含模型 一阶系统的简化 我们将研究的规范性问题是与以下形式的一阶常微分方程组： 112 3 11 221231 3 3123 1 1 112 3 1 123 1 (, , , ,) (, , , ,) (, , , ,) (, , , ,) (, , , ,) NN NN NN N NN N N fvvv v v tv vfvvvvt v fvvv v v t v f vvv v v t v f vvv v v t ? ? ? ? ?? ? ? ??? ? ??? ? ??? ? ??? = ? ??? ? ??? ? ??? ? ??? ??? ? ??? ?  "  "  " # #  "  " 或矢量形式： (,)vfvt= G G G  i v的初始条件为在t=0时的值， i v 1iN= →，许多物理系统一开始就是这个形式， 但还有一些不然。一个典型的反例是二阶常微分方程，例如： 3 ?(0) () (0) 0 vvvv v ft v =++= ? ? = ?    方程初始条件 二阶非线性常微分方程 ， ① 把这个问题转化为一个一阶系统可以简单地通过将高阶导数引进一个额外状态 量来解决。例如： 定义： 1 2 vv vv = =  则，问题①可表为： 3 221 ()vvv ft++=，其中， 1 2 ?(0) (0) 0 vv v = = 或者，重新排列为： 21 3 2 21 () vv v vvft ???? = ???? ??+ ????   112 212 (, ,) (, ,) f vvt f vvt ← ← （这是规范形式） 总之，就是把高阶系统转化为一阶的： ① 除了最高阶导数，为所有的高阶导数引入状态量 ② 用状态量代替导数，再重新排列。 线性常微分方程组的性质 某些时候，我们对线性常微分方程组系统感兴趣，举例来说： * 稳定性问题经常通过限制一个系统的均值或其它一些合理的条件来解决。 * 解析方法和数值方法比较起来，解析方法做起来更加容易一些。 * 通过线性系统可以加深对很多物理问题的理解。 我们考虑一个线性形式的规范问题： 111 222 111 () () () () NNN NNN vvft vvft A vvf vvft ??? ????? ????? ????? ????? =+ ????? ????? ????? ?????   ###   t ? ? ? ? ? ? ? ? （关于的隐含系统） i v 用矢量形式表示： ()vAvft=+ K KK  ，其中， 0 (0)vv= K K 力的这一项导致特解。考虑齐次问题: vAv= K  K ，其中， 0 (0)vv= K K 可以通过齐次问题求出特征值。 特别地，假设简并的特征值中没有重复的，那么其解就是下面的形式： 1 () j N t jj j vt e r λ α = = ∑ K K 其中： j α =常数 j λ =A的第j个特征值 j r K = A的第j个特征向量 j α， j λ，，这些都可以是复数！ j r K 我们可以把特征值分为实数部分和虚数部分： () ()ri jj j iλ λλ=+ 并且注意到： () () () () () [cos sin ] ri jjj r j ttit t ii jj eee et λλλ λ λλ = =+t ↑ 3/ 振幅由 ()r j λ控制 振荡频率由 ()i j λ控制 常微分方程组中的刚性 对于描述方程系统，刚性是一个一般的（有些模糊的）术语，它展示了宽范 围变化尺度的现象。对于我们当前感兴趣常的微分方程组的情况，也就是宽范围 变化时间尺度。 测量一个系统的刚性的一种（非常普通的）方法是把它线性化，然后计算最 大特征值的模与最小特征值的模的比，这个比值被称为谱条件数（SCN）： 谱条件数 1 1 max min j jN j jN λ λ =→ =→ ≡ 典型地，如果这个比率>1000，这个问题就被认为是刚性的。假设我们有下面的 系统： 11 22 11 0 1000 vv? ?? ???? = ???? ?? ? ???? ??   这个矩阵的特征值是-1和-1000， 刚性系统。 1000!SCN?=→ 那又怎么样呢？问题是这个系统一部分以 t e ? 缓慢衰减，另一部分以迅速衰 减。就像我们已经看到的（并且将会在接下来的一些课程中更加详细的学习 ）， 数值的稳定性通常被最快的模式（即最少的时间尺度）控制。但是我们对长时间 1000 e ? 行为更加感兴趣。所以，这导致了： max λ →控制稳定性 min λ →设置长时间行为的时间标 从我们前面在第一课中向前欧拉的例子中可知，数值运算法则的稳定性要求： max 2tλ Δ< （假设 max λ <0） max 21 500 t λ ?Δ < = ?即缓慢时标（即 min λ）的每个特征时间段500个时间步长。这样的话效 率可能会非常低！ 力项()f t K 同样能够产生长或短的时标。例如，考虑下面的问题： 例子：力的作用下的刚性问题 1000 sinvv+= t，其中，(0) 1v = 我们看下面的齐次问题： 1000 0 nn vv+ = 解得： 1000 () t n vt ae ? = ? λ =-1000 ?与外力比较“迅速”地衰减 * 在航空应用上这是非常普遍的 * 例如：体系中的声速与外力的比较由飞行器的运动决定。 所以，上面的问题也是刚性问题！对于依赖于时间的问题（这也与非线性问题有 关），一个更加普通的测量刚度的式子是： 刚性=最长时标/最短时标 向前欧拉问题 1000 sinvv tω+ = (0) 1 v = 向前欧拉问题 ：图像放大至t→0 让我们看一下向前欧拉方法是如何解决刚性力的问题的。（见下一页的附图） 结果清楚地显示了同前面相同的稳定性趋势。由2tλΔ ≤，我们得出稳定性 解（相应于Δt≤0.002），但是我们同样可以由1tλ 2< ≤+（相应于 ）得出震荡解。 0.001 0.002t<Δ ≤ 发生什么情况了？能够通过我们的直觉找到关于这个问题的一个可能 解吗？让我们回到一般的问题的标量方程： (,) dv f vt dt =，其中， 0 (0)vv= 向前欧拉问题是： 1 (,) nn nn vv f vt t + ? = + 这是一个外推法，它由时 n tt= dv dt 外推得到 1n v + ： 1 n n nn t tt dv vvt dt + = =+Δ  在时的泰勒级数，取前两项 1 (,) nn n vvtfvt + ?=+Δ n + t t n+1 t n v v n v exact (t) v n+1 有人会怀疑当Δt增加时，这样外推会有问题。 一个不同的方法是由t向回内推。例如，把v以v 1n n 1n+ 进行泰勒展开： 1 2 12 2 1 2 nn nn tt dv d v vv t t dt dt + + =?Δ Δ 去掉二次项，得： 1 1 n nn t dv vvt dt + + =+Δ 11 (,) nn nn vvtfvt Euler ++ =+Δ 向后 t+= 让我们先来看一下它是如何操作的： vv 1000 sin (0) 1v = 回想：向前欧拉方法在Δt=0.005时仍是不稳定的。 2 0 () T Evv ρ =? ∫ dt 让我们从更加近的距离看一下向后欧拉方法的稳定性。简化为一个简单的问 题： vvλ= 1 1 nn n VV V t λ + + ? ?= + 1 1 1 nn VV tλ + ?= ? + 0 1 () 1 nn VV tλ ?= ? + 如果0λ <（这自然是稳定的），那么由 1 0 1 tλ 1< ≤ ? + 推出：无论Δt多大， 的振幅是不会增长的。 n v * 这意味着Δt可以根据所需要的精度进行选择，而不是数值的稳定性。 棒极了！但是这样做的代价是什么呢？世上没有免费的午餐… * 额外的代价是每一步迭代需要更多的计算。