刚性微分方程组隐式龙格库塔方法

毕业设计

题目：刚性系统的隐式RK方法

学院：数理学院

专业名称：信息与计算科学

学号： ************

学生姓名：**

指导教师：***

2016年05月15日

摘要

本文主要介绍单步隐式Runge – Kutta方法,简要的介绍了Gauss型隐式Runge – Kutta方法、Radau型隐式Runge – Kutta方法和Lobatto型隐式Runge – Kutta方法。并利用这些基本的隐式Runge – Kutta方法来对刚性方程组进行数值求解，并将隐式Runge – Kutta方法与显式经典Runge – Kutta方法求解的结果进行对比，说明两种数值解法的优缺点。

关键词：刚性系统隐式Runge – Kutta方法单步方法Newton迭代法

Abstract

This paper mainly introduces the Implicit Runge-Kutta Methods and a simple description of Gauss implicit Runge-Kutta method ,Radau implicit Runge-Kutta method and Lobatto implicit Runge-Kutta method . These basic Implicit Runge-Kutta methods are used to solve the stiff equations. These implicit Runge-Kutta methods iare compared with the classical explicit Runge-Kutta method. This paper explain the advantages and disadvantages of the two kind of numerical methods.

Keywords: Stiff system Implicit Runge-Kutta method One step method Newton iterative method

目录

1.绪论 (1)

1.1刚性方程 (1)

1.2隐式RK方法的研究意义 (2)

1.3 RK方法的研究现状 (3)

2.单步RK方法的收敛性和稳定性 (5)

2.1单步RK方法的收敛性 (5)

2.2单步RK方法的稳定性 (6)

3.三类隐式RK方法 (8)

3.1引言 (8)

3.2 Gauss型隐式RK方法 (9)

3.3 Radau型隐式RK方法 (10)

3.4 Lobatto型隐式RK方法 (11)

4隐式RK方法的实现 (13)

4.1非线性系统的改进 (13)

4.2简化的Newton迭代法 (13)

5数值实验与结果分析 (15)

参考文献 (18)

附录 (21)

1.绪论1.1刚性方程

对于一般的线性常系数系统

y′=Ay+φ(t) A为m×m的矩阵，特征值为λi(i=1,2,⋯,m)。

定义1[23]若一个系统满足

（1）Reλi<0, i=1,2,⋯,m

（2）max

i |Reλi|min

i

|Reλi|

⁄=R≫1

其中R为刚性比，则这个系统称为刚性系统。

定义2[27]若线性系统

y′=Ay x∈[0,T]

或非线性系统

y′=f(x,y) x∈[0,T]

的矩阵A或Jacobi矩阵ðfðy

⁄的特征值λi满足

max

1≤i≤m

|Reλi|≫1

则其是刚性的。

定理1（解的存在性与唯一性）

（1）对于所有(x,y)∈D，函数f(x,y)是连续的；

（2）对于任何(x,y),(x,y∗)∈D，存在常数L，是函数满足

‖f(x,y)−f(x,y∗)‖≤L‖y−y∗‖则初值问题

{ẏ=f(x,y)(a≤x≤b) y(a)=η

有唯一解。

其中y=(y1,y2,⋯,y m)T，D={(x,y)|a≤x≤b,−∞

定义3 如果一个常微分系统的Lipscℎitz常数L很大(大于20)，则它是刚性的。

1.2隐式RK方法的研究意义

在常微分方程及常微分方程组的数值解法中，Runge – Kutta方法是目前应用最为广泛的数值解法之一，同时又具有误差小，精度高的特点。尽管显式Runge – Kutta方法能够非常准确、快速的给出大部分常微分方程组的数值解。但是在化学、自动控制电力系统等领域中，会出现一些病态的常微分方程组，也就是刚性方程组。刚性方程组对于数值解法的稳定性要求苛刻，比如方程组

{y1=−0.01y1−99.99y2, y1(0)=2

y2=−100y2, y2(0)=1将其表示为矩阵形式：

[y1

y2]=[

−0.01−99.99

0−100

]⌈

y1

y2⌉

令

A=[−0.01−99.99

0−100

]

发现A特征值为：λ1=−100，λ2=−0.01，刚性比s=|λ1|

|λ2|

=10000≫1。方程组的解为：

{y1(x)=e−100x+e−0.01x y2(x)=e−100x

解由快瞬态和慢瞬态两部分构成。

由于慢瞬态的部分，y1(x)衰减变得十分缓慢。当自变量变到x=391时，函数值还未下降到初值的1%，求解区间至少取为(0,391)。另一方面，由于快瞬态的部分，y2(x)衰减的非常快，因此步长要取得非常小。从绝对稳定性的方面来看，如果用四阶显式经典RK方法求解，绝对稳定区间要求λℎϵ(−2.78,0)，则要求ℎ< 0.0278。这样，在(0,391)上就要计算14 065步，计算量巨大，因此计算区间(0,391)内的解时，舍入误差积累会特别严重。例如取求解区间为[0,1]，用不同步长ℎ来计算y1(1)和y2(1)的值。利用四阶显式经典RK方法求解如下：

ℎy1(1)y2(1)

0.04 2.9802322e+17 2.9802322e+17

0.029.9004983e-01 1.3929556e-24

快瞬态慢瞬态