16901讲义笔记2002211 主题:*简化为一阶系统 *线性方程组系统的特性 刚度 隐含模型与显含模型 阶系统的简化 我们将研究的规范性问题是与以下形式的一阶常微分方程组 i11「f(n,n2,n…x-,nx,) f(v,V2, v3.VN-I,VN, t) fN-I(V,V2, V3. VN-,WN, t) fN(v,v2,v, 或矢量形式 =f(,1) ν的初始条件为在r=0时的v值,i=1→N,许多物理系统一开始就是这个形式, 但还有一些不然。一个典型的反例是二阶常微分方程,例如: i+y3+v=f(1) v(0)= 方程初始条件, 阶非线性常微分方程v(O)=0 把这个问题转化为一个一阶系统可以简单地通过将高阶导数引进一个额外状态 量来解决。例如: 定义 则,问题①可表为: =f(),其中 (0)=i v2(0)=0
16.901 讲义笔记 2002.2.11 主题： * 简化为一阶系统 * 线性方程组系统的特性 * 刚度 * 隐含模型与显含模型 一阶系统的简化 我们将研究的规范性问题是与以下形式的一阶常微分方程组： 1 112 3 1 2 2123 1 3 3123 1 1 112 3 1 123 1 ( , , , ,) ( , , , ,) ( , , , ,) ( , , , ,) ( , , , ,) N N N N N N N N N N N N N v fvvv v v t v f vvv v v t v f vvv v v t v N f vvv v v t v f vvv v v t − − − − − − − ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎢⎣ ⎥⎦ " " " # # " " 或矢量形式： v f vt = (,) G G G i v 的初始条件为在 t=0 时的vi 值，i =1→ N ，许多物理系统一开始就是这个形式， 但还有一些不然。一个典型的反例是二阶常微分方程，例如： 3 ( ) (0) ˆ (0) 0 v v v ft v v v + += = ⎫ ⎬ = ⎭ 方程初始条件 二阶非线性常微分方程 ， ① 把这个问题转化为一个一阶系统可以简单地通过将高阶导数引进一个额外状态 量来解决。例如： 定义： 1 2 v v v v = = 则，问题①可表为： 3 221 v v v ft ++= ( )，其中， 1 2 (0) ˆ (0) 0 v v v = =
或者,重新排列为: fv, v2, 1) 这是规范形式 2-1+f()←f(n,"2,1) 总之,就是把高阶系统转化为一阶的: ①除了最高阶导数,为所有的高阶导数引入状态量 ②用状态量代替导数,再重新排列 线性常微分方程组的性质 某些时候,我们对线性常微分方程组系统感兴趣,举例来说: *稳定性问题经常通过限制一个系统的均值或其它一些合理的条件来解决 *解析方法和数值方法比较起来,解析方法做起来更加容易一些。 通过线性系统可以加深对很多物理问题的理解。 我们考虑一个线性形式的规范问题: f(o f() 关于v的隐含系统) NN-(t) f() 用矢量形式表示 v=A+f(1),其中,v(0)=5 力的这一项导致特解。考虑齐次问题: A,其中,诉(0)= 可以通过齐次问题求出特征值。 特别地,假设简并的特征值中没有重复的,那么其解就是下面的形式: a.e 其中 =常数 2=A的第j个特征值 F=A的第j个特征向量
或者，重新排列为： 1 2 3 2 2 1 ( ) v v v v v ft ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ = ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ −−+ 112 212 ( , ,) ( , ,) f vvt f vvt ← ← （这是规范形式） 总之，就是把高阶系统转化为一阶的： ① 除了最高阶导数，为所有的高阶导数引入状态量 ② 用状态量代替导数，再重新排列。 线性常微分方程组的性质 某些时候，我们对线性常微分方程组系统感兴趣，举例来说： * 稳定性问题经常通过限制一个系统的均值或其它一些合理的条件来解决。 * 解析方法和数值方法比较起来，解析方法做起来更加容易一些。 * 通过线性系统可以加深对很多物理问题的理解。 我们考虑一个线性形式的规范问题： 1 11 2 22 1 11 ( ) ( ) ( ) ( ) N NN N NN v v ft v v ft A v vf v v ft − −− ⎡ ⎤ ⎡ ⎤⎡ ⎢ ⎥ ⎢ ⎥⎢ ⎢ ⎥ ⎢ ⎥⎢ = + ⎢ ⎥ ⎢ ⎥⎢ ⎢ ⎥ ⎢ ⎥⎢ ⎣ ⎦ ⎣ ⎦⎣ # ## t ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ （关于 的隐含系统） i v 用矢量形式表示： v Av f t = + ( ) K K K ，其中， 0 v v (0) = K K 力的这一项导致特解。考虑齐次问题: v Av =K K ，其中， 0 v v (0) = K K 可以通过齐次问题求出特征值。 特别地，假设简并的特征值中没有重复的，那么其解就是下面的形式： 1 ( ) j N t j j j vt e r λ α = = ∑ K K 其中： α j =常数 λ j =A 的第 j 个特征值 j r K = A 的第 j 个特征向量
a,,F,这些都可以是复数! 我们可以把特征值分为实数部分和虚数部分 ,=4y+i20 并且注意到 =e[cos a@t+isin aot] 个 振幅由λη)控制振荡频率由A0控制 常微分方程组中的刚性 对于描述方程系统,刚性是一个一般的(有些模糊的)术语,它展示了宽范 围变化尺度的现象。对于我们当前感兴趣常的微分方程组的情况,也就是宽范围 变化时间尺度。 测量一个系统的刚性的一种(非常普通的)方法是把它线性化,然后计算最 大特征值的模与最小特征值的模的比,这个比值被称为谱条件数(SCN): max 谱条件数≡N 典型地,如果这个比率>1000,这个问题就被认为是刚性的。假设我们有下面的 系统: 0-1000儿v2 这个矩阵的特征值是-1和-1000,→SCN=1000!→>刚性系统 那又怎么样呢?问题是这个系统一部分以e缓慢衰减,另一部分以e-0迅速衰 减。就像我们己经看到的(并且将会在接下来的一些课程中更加详细的学习) 数值的稳定性通常被最快的模式(即最少的时间尺度)控制。但是我们对长时间
α j ，λ j ， ，这些都可以是复数！ j r K 我们可以把特征值分为实数部分和虚数部分： ( ) () r i jj j λ = + λ λi 并且注意到： ( ) () ( ) ( ) ( ) [cos sin ] r i j jj r j t ti t t i i j j e ee e ti λ λλ λ λ λ = = + t ↑ 3 / 振幅由 ( )r λ j 控制 振荡频率由 ( )i λ j 控制 常微分方程组中的刚性 对于描述方程系统，刚性是一个一般的（有些模糊的）术语，它展示了宽范 围变化尺度的现象。对于我们当前感兴趣常的微分方程组的情况，也就是宽范围 变化时间尺度。 测量一个系统的刚性的一种（非常普通的）方法是把它线性化，然后计算最 大特征值的模与最小特征值的模的比，这个比值被称为谱条件数（SCN）： 谱条件数 1 1 max min j j N j j N λ λ = → = → ≡ 典型地，如果这个比率>1000，这个问题就被认为是刚性的。假设我们有下面的 系统： 1 1 2 2 1 1 0 1000 v v v v ⎡ ⎤ ⎛ ⎞ − ⎡ ⎤ ⎢ ⎥ = ⎜ ⎟⎢ ⎥ ⎣ ⎦ ⎝ ⎠ − ⎣ ⎦ 这个矩阵的特征值是-1 和-1000，⇒ = SCN 1000!→ 刚性系统。 那又怎么样呢？问题是这个系统一部分以 t e− 缓慢衰减，另一部分以 迅速衰 减。就像我们已经看到的（并且将会在接下来的一些课程中更加详细的学习 ）， 数值的稳定性通常被最快的模式（即最少的时间尺度）控制。但是我们对长时间 1000 e−
行为更加感兴趣。所以,这导致了 λ→控制稳定性 -m→设置长时间行为的时间标 从我们前面在第一课中向前欧拉的例子中可知,数值运算法则的稳定性要求 AmM<2(假设m0) m500 →即缓慢时标(即-)的每个特征时间段500个时间步长。这样的话效 率可能会非常低! 力项f(1)同样能够产生长或短的时标。例如,考虑下面的问题 例子:力的作用下的刚性问题 +1000y=sint,其中,v(0)= 我们看下面的齐次问题 i+1000n=0 解得:vn() =1000 与外力比较“迅速”地衰减 在航空应用上这是非常普遍的 *例如:体系中的声速与外力的比较由飞行器的运动决定 所以,上面的问题也是刚性问题!对于依赖于时间的问题(这也与非线性问题有 关),一个更加普通的测量刚度的式子是: 刚性=最长时标最短时标 向前欧拉问题v+1000= sin ot y(0)=1
行为更加感兴趣。所以，这导致了： λmax → 控制稳定性 λmin → 设置长时间行为的时间标 从我们前面在第一课中向前欧拉的例子中可知，数值运算法则的稳定性要求： max λ Δ < t 2 （假设λmax <0） max 2 1 500 t λ ⇒Δ < = ⇒即缓慢时标（即 λmin ）的每个特征时间段 500 个时间步长。这样的话效 率可能会非常低！ 力项 f ( )t K 同样能够产生长或短的时标。例如，考虑下面的问题： 例子：力的作用下的刚性问题 v v + = 1000 sin t ，其中，v(0) 1 = 我们看下面的齐次问题： 1000 0 n n v v + = 解得： 1000 ( ) t n v t ae− = ⇒ λ =-1000 ⇒与外力比较“迅速”地衰减 * 在航空应用上这是非常普遍的 * 例如：体系中的声速与外力的比较由飞行器的运动决定。 所以，上面的问题也是刚性问题！对于依赖于时间的问题（这也与非线性问题有 关），一个更加普通的测量刚度的式子是： 刚性=最长时标/最短时标 向前欧拉问题 v v +1000 sin = ωt v(0) 1 =
dt=0.001 0.5 dt=0019 dt=0.002 2000 向前欧拉问题:图像放大至t0
向前欧拉问题 ：图像放大至 t→0
dt=0001 dt=9019 0.05 0.05 让我们看一下向前欧拉方法是如何解决刚性力的问题的。(见下一页的附图) 结果清楚地显示了同前面相同的稳定性趋势。由|λΔl≤2,我们得出稳定性 解(相应于△01002),但是我们同样可以由1<1a2(相应于 0001<Mt≤002)得出震荡解。 发生什么情况了?能够通过我们的直觉找到关于这个问题的一个可能 解吗?让我们回到一般的问题的标量方程: f(v,r),其中,v(0) dt 向前欧拉问题是: f(v",t") 这是一个外推法,它由t=时一外推得到v+ dt
让我们看一下向前欧拉方法是如何解决刚性力的问题的。（见下一页的附图） 结果清楚地显示了同前面相同的稳定性趋势。由 λΔt ≤ 2，我们得出稳定性 解（相应于 Δt≤0.002 ），但是我们同样可以由 1< λ+t ≤ 2 （相应于 0.001 0.002 <Δ ≤t ）得出震荡解。 发生什么情况了？能够通过我们的直觉找到关于这个问题的一个可能 解吗？让我们回到一般的问题的标量方程： (,) dv f v t dt = ，其中， 0 v v (0) = 向前欧拉问题是： 1 (,) n n v v n n f v t t + − = + 这是一个外推法，它由 时n t t = dv dt 外推得到 n 1 v + ：
dv vn+△t 在t="时的泰勒级数,取前两项 v"+△f(v",t") texas o) 有人会怀疑当Mt增加时,这样外推会有问题。 一个不同的方法是由t向回内推。例如,把γ以y进行泰勒展开: 去掉二次项,得: yn+△ y+1=y+△f(vy+,rm 向后Eler 让我们先来看一下它是如何操作的 v+1000v=sntv(0)=1
1 n n n n t t t dv v vt dt + = = +Δ 在 时的泰勒级数，取前两项 1 (,) nn n v v tf v t ⇒ = +Δ + n + t t n+1 t n v v n vexact(t) v n+1 有人会怀疑当Δt 增加时，这样外推会有问题。 一个不同的方法是由t 向回内推。例如，把v 以v n 1 n n+1 进行泰勒展开： 1 2 1 2 2 1 2 n n n n t t dv d v vv t t dt dt + + = −Δ + Δ 去掉二次项，得： 1 1 n n n t dv vvt dt + + = +Δ 1 1 (,) n n nn v v tf v t Euler + + = +Δ 向后 +1 + = t 让我们先来看一下它是如何操作的： v v 1000 sin v(0) 1 =
Backward Euler:dt=0.00050.001,0.02,00021 tme 0015 0025 Backward Euler: dt=0.005 VS, 0.6 回想:向前欧拉方法在△=0005时仍是不稳定的
回想：向前欧拉方法在Δt=0.005 时仍是不稳定的
E=C(,-v)dt Farward Eule Backward Euler 让我们从更加近的距离看一下向后欧拉方法的稳定性。简化为一个简单的问 题 1-A△t 1-A△t 如果<0(这自然是稳定的,那么f1≤1推出:无论M多大,v 的振幅是不会增长的。 *这意味着M可以根据所需要的精度进行选择,而不是数值的稳定性。 棒极了!但是这样做的代价是什么呢?世上没有免费的午餐. *额外的代价是每一步迭代需要更多的计算
2 0 ( ) T E vv = − ∫ ρ dt 让我们从更加近的距离看一下向后欧拉方法的稳定性。简化为一个简单的问 题： v v = λ 1 1 n n V V n V t λ + ⇒ = − + + 1 1 1 n n V V λ t ⇒ = + − + 1 0 ( ) 1 n n V V λ t ⇒ = − + 如果λ < 0（这自然是稳定的），那么由 1 0 1 λ t < ≤1 − + 推出：无论Δt 多大， 的振幅是不会增长的。 n v * 这意味着Δt 可以根据所需要的精度进行选择，而不是数值的稳定性。 棒极了！但是这样做的代价是什么呢？世上没有免费的午餐… * 额外的代价是每一步迭代需要更多的计算