7刚性问题

预估校正法
__
• 预估
y k +1 = 2 y k − y k −1
_ 2 2 y k +1 = ( I − hA ) −1 ( hf ( xk +1 , y k +1 ) − 3 3 • 校正 __ 4 1 A y k +1 ) + y k − y k −1 ) 3 3 (7-4) __ ∂f A = ( ) ( x k +1 , y k +1 ) • 其中 ∂y 对于常系数线性方程组（7-2）上式可以简 化为 2 1 2 −1 4 yk +1 = ( I − hA) ( yk − yk −1 + hb( xk +1 )) 3 3 3 3 (7-5)
的第i个特征值， 的第 个特征值，比值 个特征值 max| Re λi |/min| Re λi | 称为刚性比，也称Stiff比。方程（1）的 称为刚性比，也称 比 方程（ ） 刚性比是100/0.01=10000 刚性比是
• 这类问题解的特点是开始阶段变化很剧烈， 这类问题解的特点是开始阶段变化很剧烈， 之后进入稳定阶段。 之后进入稳定阶段。 • 从数值方法的角度，开始阶段自然应该应 从数值方法的角度， 用精度高的Runge-Kutta方法，但是进入 方法， 用精度高的 方法 稳定阶段从数值稳定性角度， 稳定阶段从数值稳定性角度，步长仍然有 限制。 限制。 • 解决问题方法：寻找稳定的差分格式！ 解决问题方法：寻找稳定的差分格式！ • 例如隐格式。 例如隐格式。
例如：
• • • • • 求下列方程组的数值解 y’1=0.04(1-y1)-(1-y2) y1+0.0001(1-y2)2 y’2= -104y’1 +3000(1-y2)2 y1 (0)=0 , y2(0)=1, 0≤x≤100 如果用经典Runge-Kutta方法计算，必须 方法计算， 如果用经典 方法计算 h≤0.0005,若用 若用h=0.001，计算结果发散。若 若用 ，计算结果发散。 取h=0.00025，需要 ×105步。 ，需要4× • 如果在 如果在[0,0.1]区间，用h=0.00025经典 区间， 区间 经典 Runge-Kutta方法计算，在区间 方法计算， 方法计算 在区间[0.1,100]用 用 h=0.2的差分格式 的差分格式(7-5) ，则整个计算只需要 的差分格式 9×102步。 ×
y B
a
• 问题归结成非线性 方程求根问题： y(x ,t2) 方程求根问题： • y( x, t)–B=0 y(x ,tk) • 通常可以采取二分 y(x ,t3) 法、插值法、牛顿 插值法、 y(x ,t1) 法。并且把这些方 法称为打靶法。 法称为打靶法。 b x y’(a)=Байду номын сангаасk 称为试射 速度。 速度。
• 如果用经典的 如果用经典的Runge-Kutta方法求解，该方 方法求解， 方法求解 法的稳定区间是(-2.78，0)所以由第二个方 法的稳定区间是 ， 所以由第二个方 程看到，步长h应该满足 程看到，步长 应该满足 • 100h<2.78 h<0.0278 • 这样，在求解区间（0，391）上必须计算 这样，在求解区间（ ， ） 14065步，计算量很大。 步 计算量很大。
• 对于一般问题
•
y’ = Ay +b( x) (7-2) • y(x0)= y0 • 其中A是m阶方阵。 阶方阵。 • 如果系数矩阵A的特征值λi具有如下特 的特征值 征
Re λi <0 max| Re λi |/min| Re λi |是很大的数， 是很大的数， 是很大的数 则称之为刚性方程组。 则称之为刚性方程组。
9有限差分方法
• 有限差分方法的中心思想是 分化区间[a, b]， 有限差分方法的中心思想是:分化区间 分化区间 ， 例如h=(b-a)/n, xk=a+kh,k=0,1,…,n 例如 • 用差商逼近导数，离散方程和边界条件，建 用差商逼近导数，离散方程和边界条件， 立关于解函数在离散点处函数值的代数方程 组进行求解。 组进行求解。 • 例如 • y’’ = f ( x, y, y’) x∈( a, b) (8-1) ∈ • y(a)=A , y(b)=B (8-2) • 用差商代替微商离散化
8二阶边界值问题
• • • • • • • • 常微分方程二阶边界值问题一 般可写成 y’’ = f ( x, y, y’) x∈( a, b) (8-1) ∈ 并且结合下列三种边界条件之一 y(a)=A , y(b)=B (8-2) y’(a)=A , y’(b)=B (8-3) y’(a)-A0 y(a)=A1 , y’(b)+B0 y(b)=B1 (8-4) (8-4)中A0 ≥0，A0+B0 >0,它们分 中 ， 它们分 别称为第一、 三边界问题。 别称为第一、二、三边界问题。
7刚性问题 刚性问题
• 化学反应、自动控制、电力系统等领域经 化学反应、自动控制、 常遇到一类病态方程组称为刚性方程组或 方程组。 者Stiff方程组。例如 方程组 • • •
• • •
y’1 = -0.01y1 - 99.99y2 (7-1) y’2 = -100y2 y1 (0)=2 , y2(0)=1 其解为y 其解为 1 = e -0.01x + e -100x , y1= e -100x 由于e 变化很慢， ≤391 ≤391时 由于 -0.01x 变化很慢，当x≤391时， y1 ≥ y1(0)/100，即y1变化不到 变化不到1%。 ， 。
y ( xi +1 ) − 2 y ( xi ) + y ( xi −1 ) = 2 h 2 y ( xi +1 ) − y ( xi −1 ) h f [ xi , y ( xi ), − y '''(ηi )] 2h 6 2 h (4) + y (ξi ) 12
略去高阶无穷小， 略去高阶无穷小，得到相应的非线性方程 组
非线性方程组
• y’ = f ( x, y) (7-3) • y(x0)= y0 • 也存在同样的问题，只是λi表示矩阵 也存在同样的问题， ∂ f1 ∂ f1 ⋯ ∂y ∂ym 1 ∂f = ⋮ ⋮ ∂y ∂fm ∂fm ⋯ ∂ y1 ∂ym
y(xi+1)−2y(xi )+ y(xi−1) y(xi+1)−y(xi−1) + f (xi , y(xi ), ) =0 2 h 2h y0 = A yn =B,(i =1 ⋯n−1) , ,2, ,
8-1打靶法 打靶法
• 考虑第一边界问题 • y’’ = f ( x, y, y’) x∈( a, b) (8-1) ∈ • y(a)=A , y(b)=B (8-2) • 通常采取逐次逼近的方法，基本过程如下： 通常采取逐次逼近的方法，基本过程如下： 取初值问题 • y’’ = f ( x, y, y’) x∈( a, b) ∈ • y(a)=A , y’(a)=tk (8-3) • 其中 k是一个数列，假设解是 其中t 是一个数列，假设解是y(x ,tk)，而且 ， 要求y(x ,tk)→y(b)=B (k→∞ 。 →∞)。 要求 → →∞ • 这时所得的解 这时所得的解y(x)就是问题的解。如图 就是问题的解。 就是问题的解