第六章 微分方程数值解法

第6章：偏微分方程数值解法6.1对流方程【6.1.1】考虑边值问题, 01,0(0,)0,(1,)1(,0)t x x u au x t u t u t u x x=<<>ìï==íï=î如果取：2/7x D =，(0.5),1,2,3j x j x j =-D =，8/49t D =，k t k t=D 求出111123,,u u u 【解】采用Crank-Nicolson 方法()11111111211222k k k k k k k k j j j j j j j j u u u u u u u u t x ++++-+-+éù-=-++-+ëûD D 11111113k k k k k kj j j j j j u u u u u u +++-+-+-+-=-+由边界条件：(0,)0x u t =，取100k ku u x-=D ，10,0,1,k ku u k ==L (1,)1u t =，41ku =-1 1 0 0 - (1+2s) -s 0 0 -s (1+2s) -s 0 -s (1+2s) -s 0 s L L L L 101210 0 0 0 (1-2s) s 0 0 s (1-2s) s 0 s ( 1 k n n u u s u u u +-éùéùêúêúêúêúêúêú=êúêúêúêúêúêúêúêúêúëûëûL L L L L 01211-2s) s 0 1 1kn u u u u -éùéùêúêúêúêúêúêúêúêúêúêúêúêúêúêúêúëûëûL 由初始条件：021(72j j u x j ==-，1,2,3j =，212()t s x D ==D -1 1 0 0 0-1 3 -1 0 0 0 -1 3 -1 0 -1 3 -1 0 1012340 0 0 0 01 -1 1 0 00 1 -1 1 0 1 -1 1 1 u u u u u éùéùêúêúêúêúêúêú=êúêúêúêúêúêúëûëû00123 0 1 1u u u u éùéùêúêúêúêúêúêúêúêúêúêúêúêúëûëû000117u u ==，0237u =，0357u =1112327u u -=，111000123123337u u u u u u -+-=-+=，11100234235317u u u u u -+-=-+=114591u =125191u =，136991u =6.2抛物形方程【6.2.1】分别用下面方法求定解问题22(,0)4(1)(0,)(1,)0u u t x u x x x u t u t 춶=ﶶïï=-íï==ïïî01,0x t <<>（1）取0.2x D =，1/6l =用显式格式计算1i u ；(2)取0.2,0.01x t D =D =用隐式格式计算两个时间步。

第6章 常微分方程数值解法 讨论一阶常微分方程初值问题(,),,()dyf x y a x bdx y a η⎧=≤≤⎪⎪⎨⎪=⎪⎩ (6.1.1)的数值解法.数值解法可区分为两大类：(1) 单步法：此类方法在计算1n x + 上的近似值1y n + 时只用到了前一点n x 上的信息.如Euler 法，Runge-Kutta 法，Taylor 级数法就是这类方法的典型代表.(2) 多步法：此类方法在计算1yn +时，除了需要n x 点的信息外，还需要12,,n n x x -- ，等前面若干个点上的信息.线性多步法是这类方法的典型代表.离散化方法1. Taylor(台劳)展开方法2. 化导数为差商的方法3. 数值积分方法一、线性多步法基本思想：是利用前面若干个节点上()y x 及其一阶导数的近似值的线性组合来逼近下一个节点上()y x 的值. 1．一般公式的形式101',,1,,ppn in ii n i i i y a yh b y n p p +--==-=+=+∑∑其中i a ,i b 为待定常数，p 为非负整数.说明：(1)在某些特殊情形中允许任何i a 或i b 为零，但恒假设p a 和p b 不能同时全为零，此时称为1p +步法，它需要1p +个初始值01,,,.p y y y 当0p =时，定义了一类1步法，即称单步法.(2) 若10b -=，此时公式的右端都是已知的，能够直接计算出1n y +，故此时称为显式方法；若10b -≠，则公式的右端含有未知项111'(,),n n n y f x y +++=此时称其为隐式方法.2．逼近准则 准确成立：101()()'(),,1,.ppn in ii n i i i y x a y xh b y x n p p +--==-=+=+∑∑【定义 6.1】 如果对任意()r y x M =，某一线性多步法准确成立，而当()y x 为某一个1r +次多项式时，线性多步法不准确成立，则称此线性多步法是r 阶的. 注：（1）方法的阶越高，逼近效果越好. （2）1p +步法的最高阶可达 22r p =+. 3．线性多步法阶与系数的关系 局部截断误差101()()'(),,1,.ppn n in ii n i i i T y x a y xh b y x n p p +--==-=--=+∑∑()01()'()(),qq n n n q n T c y x c hy x c h y x =++++其中001011011,1[()],1{1[()()2,3,.!pi i p pi i i i p pq q q i i i i c a c i a b c i a i b q q ===--==-⎧=-⎪⎪⎪=--+⎪⎨⎪⎪⎪=--+-=⎪⎩∑∑∑∑∑【定理6.1】 线性多步法是r 阶的充分必要条件是0110,0r r C C C C +====≠称1r C +为误差常数.线性多步法是相容的：满足条件010C C ==,即0011,()1pi i ppiii i a i a b===-⎧=⎪⎪⎨⎪-+=⎪⎩∑∑∑4．线性多步法的构造方法 待定系数法：r 阶方法的系数,iia b 确定，可令010,r CC C ==== 即解下面方程得到1,0()1011()(),2,3,,01p a ii p pi a b i i i i p pq q i a q i q r i i i ⎧=∑⎪⎪=⎪⎪-+=∑∑⎪⎨==-⎪⎪⎪-⎪-+-=∑∑⎪==-⎩二、线性多步法的收敛性 记1(),pp p iii r ra rρ+-==-∑1().pi p ii r b rσ-=-=∑分别称为线性多步法的第一、第二特征多项式.()r ρ以及相应的线性多步法满足根条件：若()r ρ的所有根的模均不大于1，且模为1的根是单根。

第六章 常微分方程数值解法——RK 4法、AB 4法******（学号） *****（姓名）上机题目要求见教材P307,23题。

一、算法原理题目要求采用RK 4法和AB 4法求解最简单的常微分方程初值问题(,),()y f x y a x by a η'=≤≤⎧⎨=⎩ (1)为求解式(1)，采用离散化方法，就是寻求解)(x y 在区间],[b a 上的一系列点<<<<<n x x x x 321上的近似值 ,,,,21n y y y 。

记1(1,2,)i i i h x x i -=-=表示相邻两个节点的间距，称为步长。

求微分方程数值解的主要问题：(1) 如何将微分方程(,)y f x y '=离散化,并建立求其数值解的递推公式; (2) 递推公式的局部截断误差、数值数n y 与精确解)(n x y 的误差估计； (3) 递推公式的稳定性与收敛性. a) Runge-Kutta 方法基本思想：通过在1[,]i i x x +多预报几个点求斜率，并将其加权平均作为k *的近似值，以此构造更高精度的计算公式。

如果每步计算四次函数 的值，完全类似的，可以导出局部截断误差为)(5h O 的四阶Runge-Kutta 公式(RK 4)：1123412132431(22),6(,),(,),221(,),22(,).n n n n n n n n n n y y k k k k k f x y h h k f x y k h k f x h y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ (2)b) Adams 显式公式Runge-Kutta 方法是单步法，计算1+n y 时，只用到n y ， 而已知信息1-n y 、2-n y 等没有被直接利用。

可以设想如果充分利用已知信息1-n y ，2-n y ，…来计算1+n y ，那么不但有可能提高精度，而且大大减少了计算量，这就是构造所谓线性多步法的基本思想。

第六章常微分方程的数值解法第六章常微分方程的数值解法在自然科学研究和工程技术领域中，常常会遇到常微分方程的求解问题。

传统的数学分析方法仅能给出一些简单的、常系数的、经典的线性方程的解析表达式，不能处理复杂的、变系数的、非线性方程，对于这些方面的问题，只能求诸于近似解法和数值解法。

而且在许多实际问题中，确确实实并不总是需要精确的解析解，往往只需获得近似的解或者解在若干个点上的数值即可。

在高等数学课程中介绍过的级数解法和逐步逼近法，能够给出解的近似表达式，这一类方法称为近似解法。

还有一类方法是通过计算机来求解微分方程的数值解，给出解在一些离散点上的近似值，这一类方法称作为数值方法。

本章主要介绍常微分方程初值问题的数值解法，包括Euler 方法、Runge-Kutta 方法、线性多步法以及微分方程组与高阶微分方程的数值解法。

同时，对于求解常微分方程的边值问题中比较常用的打靶法与有限差分法作了一个简单的介绍。

§1 基本概念1.1 常微分方程初值问题的一般提法常微分方程初值问题的一般提法是求解满足如下条件的函数，，b x a x y ≤≤)(=<<=α)(),(a y bx a y x f dxdy， (1.1) 其中),(y x f 是已知函数，α是给定的数值。

通常假定上面所给出的函数),(y x f 在给定的区域},),{(+∞<≤≤=yb x a y x D 上面满足如下条件：(1) 函数),(y x f 在区域D 上面连续；(2) 函数),(y x f 在区域D 上关于变量y 满足Lipschitz(李普希茨)条件：212121,),(),(y y b x a y y L y x f y x f ?≤≤?≤?，， (1.2)其中常数L 称为Lipschitz(李普希茨)常数。

由常微分方程的基本理论可以知道，假如(1.1)中的),(y x f 满足上面两个条件，则常微分方程初值问题(1.1)对于任意给定的初始值α都存在着唯一的解，，b x a x y ≤≤)(并且该唯一解在区间[a,b]上是连续可微的。

第六章 常微分方程的数值解法§6.0 引言§6.1 算法构造的主要途径 §6.2 Runge-Kutta Method 算法 §6.3 线性多步法§6.4 线性多步法的一般形式 §6.5 算法的稳定性、收敛性§6.0 引 言1． 主要考虑如下的一阶常微分方程初值问题的求解:()()00,dyf x y dx y x y ⎧=⎪⎪⎨=⎪⎪⎩微分方程的解就是求一个函数y=y(x)，使得该函数满足微分方程并且符合初值条件。

2. 例如微分方程:xy '-2y=4x ；初始条件: y(1)=-3。

于是可得一阶常微分方程的初始问题24(1)3yy x y ⎧'=+⎪⎨⎪=-⎩。

显然函数y(x)=x 2-4x 满足以上条件,因而是该初始问题的微分方程的解。

3. 但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示。

因此，只能依赖于数值方法去获得微分方程的数值解。

4． 微分方程的数值解：设微分方程问题的解y(x)的存在区间是[a,b ]，初始点x 0=a ，将[a,b ]进行划分得一系列节点x 0 , x 1 ,...,x n ，其中a= x 0< x 1<…< x n =b 。

y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点x k 的近似值y(x k )，即y≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。

如图所示:a bx 0x 1x 2...§6.1 算法构造的主要途径1 欧拉公式1.1 构造的思想：微分方程初值问题: ()()00,dyf x y dx y x y ⎧=⎪⎨⎪=⎩ 利用差商代替一阶导数，即10()()x x y x y x dy dxh =-≈，则1000()()(,)y x y x f x y h -≈。