08第二版 第八章 常微分方程数值解法
5 1.0 2.9766 3.4366 0.4600
（2）由隐式欧拉公式(8.6)，得
yi yi1 0.2(xi yi ), i 1, 2, ,5
整理，得
yi (0.2xi yi1) / 0.8, i 1, 2, ,5
计算结果见下表 (表8-2)：
i xi yi y( xi ) y(xi) yi
,
yi1
)
f (xi, yi )],
i
1, 2,
, n (8.7)
上式称为常微分方程初值问题(8.1)～(8.2)的梯形公式．
另外，将微分程(8.1)两端从 xi1 到 xi 积分，得
xi y(x)dx xi f (x, y(x))dx
xi1
xi1
上式右端积分用梯形公式近似，即
xi xi1
y(x) 在自变量 x的一系列离散节点
a x0 x1 x2 xn1 xn b
上的近似值
y0, y1, y2, , yn1, yn
这些近似值称为初值问题(8.1)～(8.2)的数值解.
相邻两节点的间距 hi xi xi1 (i 1, 2, , n) 称为步长 ，通常在计算上采用相等的步长hi h (i 1, 2, , n) ， 这时等 距节点 xi x0 ih, (i 1,2, ,n) .求解过程是顺着节点排列的顺序
8.1.3 微分方程单步法的局部截断误差与阶
前面几节介绍了求解初值问题(8.1)～(8.2)数值解的几种方
法．显然，各种数值方法得到的数值解 yi 与解析解 y(xi )之间
的差异各不相同．称
ei y(xi ) yi 为某方法在 xi 处的整体截断误差．显然，该误差依赖于前面
xi1, xi2 , , x0
由梯形公式(8.7)，得
yi yi1 0.1[(xi1 yi1) (xi yi )] i 1, 2, ,5
整理，得
yi [0.1(xi1 xi ) 1.1yi1] / 0.9 i 1, 2, ,5
计算结果如下：
i
0
xi yi y( xi ) y(xi) yi
0 1.0000 1.0000
y(x) 2ex x 1 相比较.
解 根据已知条件，有
n 5, h 0.2
x0 0, y0 1 xi ih, i 1, 2, , n
又
f (x, y) x y
（1）由显式欧拉公式(8.5)，得
yi yi1 0.2(xi1 yi1), i 1, 2, , 5
整理，得
yi 0.2xi1 1.2yi1,i 1, 2, ,5
y(x) x y, 0 x 1
y(0)
1
取 h 0.2 ，试用改进欧拉公式求其数值解，并与精确解
y(x) 2ex x 1 相比较.
解 根据已知条件，有
n 5, h 0.2
x0 0, y0 1
xi ih, i 1, 2, , n
又
f (x, y) x y
根据改进的欧拉公式(8.8)，有
一步一步的向前推进，初值问题的数值解法主要建立这种递推 公式．
求解常微分方程初值问题的数值解法可分为两类：单步法和多 步法．单步法在计算 yi 时只用到前面一步 yi1 处的信息，如
欧拉(Euler)法、龙格—库塔(Runge-Kutta)法．多步法在计算 yi 时用到前面多步 yi1, yi2 , 处的信息，常用的主要是线性多
，又可得
yi yi1 hf (xi1, yi1), i 1, 2, , n
（8.5）
及
yi yi1 hf (xi , yi ), i 1, 2, , n （8.6）
(8.5)、(8.6)式称为欧拉公式．其中(8.5)式称为显式欧拉公式， (8.6)式称为隐式欧拉公式．利用欧拉公式求解初值问题(8.1)～ (8.2)数值解的方法称为欧拉法．
第八章常微分方程数值解法
• 8.1 欧拉法 • 8.2 龙格-库塔法 • 8.3 线性多步法 • 8.4 一阶微分方程组与高阶微分方程的数值解法
在科学研究和工程技术中，常会遇到常微分方程或常微分 方程组的求解问题．我们知道，除了几种简单类型的常微分方 程外，要找出解的解析表达式是极其困难的，甚至是不可能 的．因此，研究各种类型常微分方程的数值解法是很有必要 的．本章主要考虑一阶常微分方程初值问题
步法．
单步法与多步法都有显式格式与隐式格式之分.
§8.1欧拉法
8.1.1 欧拉公式
设初值问题(8.1)～(8.2)的精确解为y(x)在区间[a,b] 上取一
系列等距节点
a x0 x1 x2 xn b
根据导数定义，可以用
y(xi ) y(xi1) h
近似代替
y(xi1)
或
y(xi （) 图8-1），其中 h
f (x, y1) f (x, y2) L y1 y2 .
对任何 y1, y2 均成立，则上述初值问题(7.1)～(7.2)的解存在、 惟一，且连续依赖于初值 y0．
由于李普希茨条件较难验证，因此在实际应用中，常用函
数 f (x, y) 在所考虑区域 G 上对 y 存在连续偏导数条件代替．事
实上，若 f (x, y) 在有解闭区域 G 上对 y 存在连续偏导数，则 f
计算结果见下表 (表8-1)：
i xi yi y(xi ) y(xi ) yi
0 0 1.0000 1.0000 0.0000
表8-1
1
2
3
0.2
0.4
0.6
1.2000 1.4800 1.8560
1.2428 1.5836 2.0442
0.04.6511 0.3039
此当函数 f (x, y) 比较复杂时从(8.7)式中解出yi 还是十分困难.
能否找到一种公式，既能提高精度，又可避免从代数方程
中求解 yi 这一繁琐过程呢？改进的欧拉公式正好满足这些要
求．
在实际应用中，首先利用显式欧拉公式(8.5)计算 yi 的 一个初值 yi (称为预测值)，再用梯形公式对yi 进行校正
b a .于是， n
根据微分方程(8.1)可得
y(xi ) y(xi1) h
f
(xi1, yi1),
（8.3）
i 1, 2, , n
或
y(xi ) y(xi1) h
f (xi , yi ),
i 1, 2, , n
（8.4）
设 yi ，yi1 分别为精确解 y(xi1)， y(xi ) 的近似值,则由(8.3)、(8.4)式
y(i1) h2 ,
2!
i1 在 xi1 与 xi 之间
将 y(x) 在 xi 点处展开，有
y(x)
y(xi )
y(xi )(x
xi )
y( ) (x
2!
xi )2 ,
在x与 xi 之间
令 x xi1 ，得
y( xi 1 )
y(xi )
y(xi )h
y(i )
2!
h2
y(xi )
f
(xi , y(xi ))h
dy
dx
f (x, y), a
xb
y(a) y0
（8.1） （8.2）
在区间[a, b]上的解. 为了保证初值问题(8.1)～(8.2)解的存在和惟 一性，给出如下定理(证明略).
定理8.1 若 f (x, y) 连续，且关于 y 满足李普希兹(Lipschitz)条
件，即存在常数 L 0 ，使
各步的计算．为了简化分析，我们仅讨论向前推进一步所产 生的误差，因此有如下定义．
定义8.1 利用数值方法求解微分方程的近似解时，假定 yi1 的 计算没有误差，即 yi1 y(xi1) ，则称 y(xi ) yi 为局部截断误
差．
定义8.2 如果对某种数值方法，存在正整数p,使其局部截断误
8.1.2 改进欧拉公式
根据(8.3),(8.4)式，令
y(xi ) y(xi1) f (xi1, yi1) f (xi , yi )
h
2
整理，得
yi
yi1
h[ 2
f
(xi1,
yi1)
f
( xi ,
yi )]
所以，通过取显示欧拉公式和隐式欧拉公式的算术平均值，
得
yi
yi1
h 2
[
f
(
xi1
y(i ) h2 ,
2!
在 xi1 与 xi 之间
当h xi xi1 充分小时，忽略二阶导数项，并以近似值 yi1, yi 代替精确值 y(xi1), y(xi )，同样可得欧拉公式(8.5)、(8.6)．
例8.1 设有初值问题
y(x) x y, 0 x 1
y(0)
1
取 h 0.2 ，试用梯形公式求其数值解，并与精确解
0
1 0.2 1.2444 1.2428 -0.0016
2 0.4 1.5877 1.5836 -0.0041
3 0.6 2.0516 2.0442 -0.0074
4 0.8 2.6630 2.6511 -0.0119
5 1.0 3.4548 3.4366 -0.0182
将例8.2与例8.1 (1)、(2)对比，可以看出梯形公式比欧拉公 式的精度明显提高．但由于梯形公式(8.7)仍然是隐式公式，因
，就可得到所谓改进的欧拉公式
yi
yi1
hf
( xi 1 ,
yi 1 )
,i 1, 2,
h
yi
yi1
[ 2
f
( xi 1 ,
yi1)
f
(xi ,
yi )]
, n （8.8）
(8.8)式中的第一式称为预测算式，用来得到精度较低的预测值 ；第二式称为校正算式，用来得到较高精度的校正值．
常微分方程的数值解法
常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。
然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。
在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。
这些方法可以给出解的近似表达式，通常称为近似解析方法。
还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。
利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题的精确解记为y (x )。
数值方法的基本思想是：在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。
这些h i 可以不相等，但一般取成相等的，这时na b h -=。
在这些节点上采用离散化方法，（通常用数值积分、微分。
泰勒展开等）将上述初值问题化成关于离散变量的相应问题。
把这个相应问题的解y n 作为y (x n )的近似值。
这样求得的y n 就是上述初值问题在节点x n 上的数值解。
一般说来，不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法 1.欧拉法1．对常微分方程初始问题(9.2))((9.1) ),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
常微分方程的数值解法的原理
常微分方程的数值解法的原理
常微分方程的数值解法是一种解决常微分方程初值问题的方法，其基本原理是：先取自变量的一系列离散点，把微分问题离散化，求出离散问题的数值解，以此作为微分问题解的近似。
具体步骤如下：
1、假定解存在且唯一，解函数y（x）及右端函数?（x，y）具有所需的光滑
程度。
2、取步长h＞0，以h剖分区间【α，b】，令xi=α+ih，把微分方程离散化成
一个差分方程。
3、以y（x）表微分方程初值问题的解，以yi表差分问题的解，就是近似解
的误差，称为全局误差。
4、设计各种离散化模型，求出近似解，估计误差以及研究数值方法的稳定性
和收敛性等。
第8章常微分方程数值解法
的解为
y ( x) e
x2
x 0
e dt
t2
但要计算它的值，还需要用数值积分的方法。如果要 对许多个 x 值计算解 y(x) 的近似值，那么工作量非常大。况 且实际计算不一定要求解析表达式，而是只需求在某些点 上满足精度的解的近似值或解的近似表达式就可以了。
由于高阶常微分方程可以转化为一阶常微分方程组，因 此，为了不失一般性，本章主要介绍一类一阶常微分方程初 值问题
的解来近似微分方程初值问题(8.2)的解，其 中 h (b- a) / 2 ，式(8.3)也称为欧拉公式。
欧拉法的几何意义是用一条自点 ( x0 , y0 ) 出发的 折线去逼近积分曲线 y f (x) ，如图8.1所示。 因此，这种方法又称为折线法。显然，欧拉法 简单地取折线的端点作为数值解，精度非常差。
float euler(float x0,float xn,float y0,int N) { float x,y,h; int i; x=x0; y=y0; h=(xn-x0)/(float)N; /* 计算步长 */ for(i=1;i<=N;i++) /* 欧拉公式 */ { y=y+h*func(x,y); x=x0+i*h; } return(y); }
8.4 龙格—库塔（Runge-Kutta）法 8.4.1 龙格—库塔法的基本思想
在欧拉法 yi 1 yi h f ( xi , yi ) (i 0,1,) 中，用解函数 y f (x) 在 点 x i 处的斜率 f ( xi , y i ) 计算从 yi 到 y i 1 的增量，y i 1 的表达式 与 y( xi 1 ) 的Taylor展开式的前二项相等，使方法只有一阶精度。 改进的欧拉法用两个点 x i ，x i 1 处的斜率 f ( xi , y i )、f ( xi 1 , yi 1 ) 的平均值计算增量，使方法具有二阶精度，即 y i 1 的表达式 与 y( xi 1 ) 的Taylor展开式的前三项相等。 由此龙格和库塔提出了一种间接地运用Taylor公式的方法， y (x) 即利用 在若干个待定点上的函数值和导数值做出线性组 合式，选取适当系数使这个组合式进行Taylor展开后与 y( xi 1 ) 的Taylor展开式有较多的项达到一致，从而得出较高阶的数 值公式，这就是龙格—库塔法的基本思想。
第八章常微分方程的数值解法
y( xn1 )
15
Euler法的收敛性
称初值问题(8.1.1)的数值解法是收敛的，如:
h0 ( n )
lim yn y ( x)
其中： x xn x0 nh , x [ x0 , b]
16
例考察以下初值问题Euler法的收敛性
dy y dx y (0)＝y0 ( 0)
★
可得： h (k ) ( k 1) y y | f ( xn 1 , yn ) f ( x , y 1 n 1 n 1 ) | 2 hL ( k ) hL k 1 (1) ( k 1) (0) | yn 1 yn 1 | ( ) | yn 1 yn 1 | 2 2 hL k 1 ( k 1) 从而 : lim( ) 0 , 故有 lim yn 1 y n 1 。 k 2 k
★
由y0＝y( x0 ), 假定yn＝y( xn ), 往证：
y0 yn 1 y ( xn 1 ) xn 1； x0
14
证明
yn yn1 yn hf ( xn , yn ) yn h xn 1 1 yn (1 h ) y( xn )(1 h ) xn xn y0 y0 1 xn (1 h ) ( xn h) x0 xn x0 y0 xn 1 x0
8
局部截断误差
假设第n步在点xn的值计算没有误差,即yn y( xn ), 由单步法计算出yn1 , 则
Tn1 y( xn1 ) yn1 称为点xn1上的局部截断误差.
从初值y( x0 ) y0出发，由单步法显式或隐式 逐步计算，得xn 1的值yn 1 , 则
n1 y( xn1 ) yn1
常微分方程的数值解法
常微分方程的数值解法1. 引言常微分方程是自变量只有一个的微分方程，广泛应用于自然科学、工程技术和社会科学等领域。
由于常微分方程的解析解不易得到或难以求得，数值解法成为解决常微分方程问题的重要手段之一。
本文将介绍几种常用的常微分方程的数值解法。
2. 欧拉方法欧拉方法是最简单的一种数值解法，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上假设解函数为线性函数，即通过给定的初始条件在每个子区间上构造切线；- 使用切线的斜率（即导数）逼近每个子区间上的解函数，并将其作为下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
3. 改进的欧拉方法改进的欧拉方法是对欧拉方法的一种改进，主要思想是利用两个切线的斜率的平均值来逼近每个子区间上的解函数。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上构造两个切线，分别通过给定的初始条件和通过欧拉方法得到的下一个初始条件；- 取两个切线的斜率的平均值，将其作为该子区间上解函数的斜率，并计算下一个子区间的初始条件；- 重复上述过程直至达到所需的精度。
4. 二阶龙格-库塔方法二阶龙格-库塔方法是一种更为精确的数值解法，其基本思想是通过近似计算解函数在每个子区间上的平均斜率。
具体步骤如下： - 将自变量的区间等分为n个子区间；- 在每个子区间上计算解函数的斜率，并以该斜率的平均值近似表示该子区间上解函数的斜率；- 利用该斜率近似值计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
5. 龙格-库塔法（四阶）龙格-库塔法是目前常用的数值解法之一，其精度较高。
四阶龙格-库塔法是其中较为常用的一种，其具体步骤如下：- 将自变量的区间等分为n个子区间；- 在每个子区间上进行多次迭代计算，得到该子区间上解函数的近似值；- 利用近似值计算每个子区间上的斜率，并以其加权平均值逼近解函数的斜率；- 计算下一个子区间的初始条件，并进一步逼近解函数；- 重复上述过程直至达到所需的精度。
常微分方程的数值解法(欧拉法、改进欧拉法、泰勒方法和龙格库塔法)
［例1］用欧拉方法与改进的欧拉方法求初值问题h 的数值解。
在区间[0,1]上取0.1［解］欧拉方法的计算公式为x0=0;y0=1;x(1)=0.1;y(1)=y0+0.1*2*x0/(3*y0^2);for n=1:9x(n+1)=0.1*(n+1);y(n+1)=y(n)+0.1*2*x(n)/(3*y(n)^2);end;xy结果为x =Columns 1 through 80.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 100.9000 1.0000y =Columns 1 through 81.0000 1.0067 1.0198 1.0391 1.0638 1.0932 1.1267 1.1634 Columns 9 through 101.2028 1.2443改进的欧拉方法其计算公式为本题的精确解为()y x=x0=0;y0=1;ya(1)=y0+0.1*2*x0/(3*y0^2);y(1)=y0+0.05*(2*x0/(3*y0^2)+2*x0/(3*ya^2));for n=1:9x(n+1)=0.1*(n+1);ya(n+1)=ya(n)+0.1*2*x(n)/(3*ya(n)^2);y(n+1)=y(n)+0.05*(2*x(n)/(3*y(n)^2)+2*x(n+1)/(3*ya(n+1)^2));end;xy结果为x =Columns 1 through 80.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 100.9000 1.0000y =Columns 1 through 81.0000 1.0099 1.0261 1.0479 1.0748 1.1059 1.1407 1.1783 Columns 9 through 101.2183 1.2600［例2］用泰勒方法解x＝0.1, 0.2, …, 1.0处的数值解，并与精确解进行比较。
常微分方程组数值解法
常微分方程组数值解法一、引言常微分方程组是数学中的一个重要分支，它在物理、工程、生物等领域都有广泛应用。
对于一些复杂的常微分方程组，往往难以通过解析方法求解，这时候数值解法就显得尤为重要。
本文将介绍常微分方程组数值解法的相关内容。
二、数值解法的基本思想1.欧拉法欧拉法是最基础的数值解法之一，它的思想是将时间连续化，将微分方程转化为差分方程。
对于一个一阶常微分方程y'=f(x,y)，其欧拉公式为：y_{n+1}=y_n+hf(x_n,y_n)其中h为步长，x_n和y_n为第n个时间点上x和y的取值。
2.改进欧拉法改进欧拉法是对欧拉法的改良，其公式如下：y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_n+hf(x_n,y_n))] 3.四阶龙格-库塔方法四阶龙格-库塔方法是目前最常用的数值解法之一。
其公式如下：k_1=f(x_n,y_n)k_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}k_1)k_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}k_2)k_4=f(x_n+h,y_n+hk_3)y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)其中，k_i为中间变量。
三、常微分方程组的数值解法1.欧拉法对于一个二阶常微分方程组：\begin{cases} y'_1=f_1(x,y_1,y_2) \\ y'_2=f_2(x,y_1,y_2)\end{cases}其欧拉公式为：\begin{cases} y_{n+1,1}=y_{n,1}+hf_1(x_n,y_{n,1},y_{n,2}) \\y_{n+1,2}=y_{n,2}+hf_2(x_n,y_{n,1},y_{n,2}) \end{cases}其中，x_n和y_{n,i}(i=1, 2)为第n个时间点上x和y_i的取值。
常微分方程的数值解法
常微分方程的数值解法常微分方程是研究变量的变化率与其当前状态之间的关系的数学分支。
它在物理、工程、经济等领域有着广泛的应用。
解常微分方程的精确解往往十分困难甚至不可得，因此数值解法在实际问题中起到了重要的作用。
本文将介绍常见的常微分方程的数值解法，并比较其优缺点。
1. 欧拉方法欧拉方法是最简单的数值解法之一。
它基于近似替代的思想，将微分方程中的导数用差商近似表示。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）根据微分方程的定义使用近似来计算下一个点的值。
欧拉方法的计算简单，但是由于误差累积，精度较低。
2. 改进欧拉方法为了提高欧拉方法的精度，改进欧拉方法应运而生。
改进欧拉方法通过使用两个点的斜率的平均值来计算下一个点的值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）根据微分方程的定义使用近似来计算下一个点的值。
改进欧拉方法相较于欧拉方法而言，精度更高。
3. 龙格-库塔法龙格-库塔法（Runge-Kutta）是常微分方程数值解法中最常用的方法之一。
它通过迭代逼近精确解，并在每一步中计算出多个斜率的加权平均值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）计算各阶导数的导数值。
（4）根据权重系数计算下一个点的值。
与欧拉方法和改进欧拉方法相比，龙格-库塔法的精度更高，但计算量也更大。
4. 亚当斯法亚当斯法（Adams）是一种多步法，它利用之前的解来近似下一个点的值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）通过隐式或显式的方式计算下一个点的值。
亚当斯法可以提高精度，并且比龙格-库塔法更加高效。
5. 多步法和多级法除了亚当斯法，还有其他的多步法和多级法可以用于解常微分方程。
多步法通过利用多个点的值来逼近解，从而提高精度。
而多级法则将步长进行分割，分别计算每个子问题的解，再进行组合得到整体解。
常微分方程的数值解法全文
第8章常微分方程的数值解法8.4单步法的收敛性与稳定性8.4.1相容性与收敛性上面所介绍的方法都是用离散化的方法，将微分方程初值问题化为差分方程初值问题求解的.这些转化是否合理？即当h →∞时，差分方程是否能无限逼近微分方程，差分方程的解n y 是否能无限逼近微分方程初值问题的准确解()n y x ，这就是相容性与收敛性问题.用单步法(8.3.14)求解初值问题(8.1.1)，即用差分方程初值问题100(,,)()n n n n y y h x y h y x y ϕ+=+⎧⎨=⎩(8.4.1)的解作为问题(8.1.1)的近似解，如果近似是合理的，则应有()()(,(),)0 (0)y x h y x x y x h h hϕ+--→→(8.4.2)其中()y x 为问题(8.1.1)的精确解.因为0()()lim ()(,)h y x h y x y x f x y h→+-'==故由(8.4.2)得lim (,,)(,)h x y h f x y ϕ→=如果增量函数(,(),)x y x h ϕ关于h 连续，则有(,,0)(,)x y f x y ϕ=(8.4.3)定义8.3如果单步法的增量函数(,,)x y h ϕ满足条件(8.4.3)，则称单步法(8.3.14)与初值问题(8.1.1)相容.通常称(8.4.3)为单步法的相容条件.满足相容条件(8.4.3)是可以用单步法求解初值问题(8.1.1)的必要条件.容易验证欧拉法和改进欧拉法均满足相容性条件.一般地，如果单步法有p 阶精度（1p ≥），则其局部截断误差为[]1()()(,(),)()p y x h y x h x y x h O h ϕ++-+=上式两端同除以h ，得()()(,,)()p y x h y x x y h O h hϕ+--=令0h →，如果(,(),)x y x h ϕ连续，则有()(,,0)0y x x y ϕ'-=所以1p ≥的单步法均与问题(8.1.1)相容.由此即得各阶龙格－库塔法与初值问题(8.1.1)相容.定义8.4一种数值方法称为是收敛的，如果对于任意初值0y 及任意固定的(,]x a b ∈，都有lim () ()n h y y x x a nh →==+其中()y x 为初值问题(8.1.1)的精确解.如果我们取消局部化假定，使用某单步法公式，从0x 出发，一步一步地推算到1n x +处的近似值1n y +.若不计各步的舍入误差，而每一步都有局部截断误差，这些局部截断误差的积累就是整体截断误差.定义8.5称111()n n n e y x y +++=-为某数值方法的整体截断误差.其中()y x 为初值问题(8.1.1)的精确解，1n y +为不计舍入误差时用某数值方法从0x 开始，逐步得到的在1n x +处的近似值（不考虑舍入误差的情况下，局部截断误差的积累）.定理8.1设单步法(8.3.14)具有p 阶精度，其增量函数(,,)x y h ϕ关于y 满足利普希茨条件，问题(8.1.1)的初值是精确的，即00()y x y =，则单步法的整体截断误差为111()()p n n n e y x y O h +++=-=证明由已知，(,,)x y h ϕ关于y 满足利普希茨条件，故存在0L >，使得对任意的12,y y 及[,]x a b ∈，00h h <≤，都有1212(,,)(,,)x y h x y h L y y ϕϕ-≤-记1()(,(),)n n n n y y x h x y x h ϕ+=+，因为单步法具有p 阶精度，故存在0M >，使得1111()p n n n R y x y Mh ++++=-≤从而有111111111()()()(,(),)(,,)()(,(),)(,,)n n n n n n n p n n n n n n p n n n n n n e y x y y x y y y Mh y x h x y x h y h x y h Mh y x y h x y x h x y h ϕϕϕϕ+++++++++=-≤-+-≤++--≤+-+-1(1)p nMh hL e +≤++反复递推得11111101110(1)(1)1(1)(1)(1)(1)1(1)p p n n n p n n p n e Mh hL Mh hL e hL hL Mh hL e hL Mh hL e hL+++-+++++⎡⎤≤++++⎣⎦⎡⎤≤+++++++⎣⎦+-≤++因为00()y x y =，即00e =，又(1)n h b a +≤-，于是ln(1)1()(1)(1)b a b a hL n L b a h h hL hL e e --++-+≤+=≤所以()11()p L b a p n M e h e O h L -+⎡⎤≤-=⎣⎦推论设单步法具有p (1p ≥)阶精度，增量函数(,,)x y h ϕ在区域G ：, , 0a x b y h h ≤≤-∞<<+∞≤≤上连续，且关于y 满足利普希茨条件，则单步法是收敛的.当(,)f x y 在区域:,D a x b y ≤≤-∞<<+∞上连续，且关于y 满足利普希茨条件时，改进欧拉法，各阶龙格－库塔法的增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，因而它们都是收敛的.关于单步法收敛的一般结果是：定理8.2设增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，则单步法收敛的充分必要条件是相容性条件(8.4.3).8.4.2稳定性稳定性与收敛性是两个不同的概念，收敛性是在假定每一步计算都准确的前提下，讨论当步长0h →时，方法的整体截断误差是否趋于零的问题.而稳定性则是讨论舍入误差的积累能否对计算结果有严重影响的问题.定义8.6若一种数值方法在节点值n y 上有一个大小为δ的扰动，于以后各节点()m y m n >上产生的偏差均不超过δ，则称该方法是稳定的.我们以欧拉法为例进行讨论.假设由于舍入误差，实际得到的不是n y 而是n n n y y δ=+，其中n δ是误差.由此再计算一步，得到1(,)n n n n y y hf x y +=+把它与不考虑舍入误差的欧拉公式相减，并记111n n n y y δ+++=-，就有[]1(,)(,)1(,)n n n n n n y n nh f x y f x y hf x δδηδ+⎡⎤=+-=+⎣⎦其中y f f y∂=∂.如果满足条件1(,)1y n hf x η+≤，(8.4.4)则从n y 到1n y +的计算，误差是不增的，可以认为计算是稳定的.如果条件(8.4.4)不满足，则每步误差将增大.当0y f >时，显然条件(8.4.4)不可能满足，我们认为问题本身具有先天的不稳定性.当0y f <时，为了满足稳定性要求(8.4.4)，有时h 要很小.一般的，稳定性与方法有关，也与步长h 的大小有关，当然也与方程中的(,)f x y 有关.为简单起见，通常只考虑数值方法用于求解模型方程的稳定性，模型方程为y y λ'=(8.4.5)其中λ为复数.一般的方程可以通过局部线性化转化为模型方程，例如在(,)x y 的邻域(,)(,)(,)()(,)()x y y f x y f x y f x y x x f x y y y '==+-+-+略去高阶项，再作变量替换就得到u u λ'=的形式.对于模型方程(8.4.5)，若Re 0λ>，类似以上分析，可以认为方程是不稳定的.所以我们只考虑Re 0λ<的情形，这时不同的数值方法可能是数值稳定的或者是数值不稳定的.当一个单步法用于试验方程y y λ'=，从n y 计算一步得到1()n n y E h y λ+=(8.4.6)其中()E h λ依赖于所选的方法.因为通过点(,)n n x y 试验方程的解曲线（它满足,()n n y y y x y λ'==）为[]exp ()n n y y x x λ=-，而一个p 阶单步法的局部截断误差在()n n y x y =时有1111()()p n n n T y x y O h ++++=-=，所以有1exp()()()p n n y h E h y O h λλ+-=(8.4.7)这样可以看出()E h λ是h e λ的一个近似值.由(8.4.6)可以看到，若n y 计算中有误差ε，则计算1n y +时将产生误差()E h λε，所以有下面定义.定义8.7如果(8.4.6)式中，()1E h λ<，则称单步法(8.3.14)是绝对稳定的.在复平面上复变量h λ满足()1E h λ<的区域，称为方法(8.3.14)的绝对稳定区域，它与实轴的交称为绝对稳定区间.在上述定义中，规定严格不等式成立，是为了和线性多步法的绝对稳定性定义一致.事实上，()1E h λ=时也可以认为误差不增长.(1)欧拉法的稳定性欧拉法用于模型方程(8.4.5)，得1(1)n n y h y λ+=+，所以有()1E h h λλ=+.所以绝对稳定条件是11h λ+<，它的绝对稳定区域是h λ复平面上以(1,0)-为中心的单位圆，见图8.3.而λ为实数时，绝对稳定区间是(2,0)-.Im()h λRe()h λ2-1-O 图8.3欧拉法的绝对稳定区域(2)梯形公式的稳定性对模型方程，梯形公式的具体表达式为11()2n n n n h y y y y λλ++=++，即11212n nh y y h λλ++=-，所以梯形公式的绝对稳定区域为12112h h λλ+<-.化简得Re()0h λ<，因此梯形公式的绝对稳定区域为h λ平面的左半平面，见图8.4.特别地，当λ为负实数时，对任意的0h >，梯形公式都是稳定的.Im()h λRe()h λO 图8.4梯形公式的绝对稳定区域(3)龙格－库塔法的稳定性与前面的讨论相仿，将龙格－库塔法用于模型方程(8.4.5)，可得二、三、四阶龙格－库塔法的绝对稳定区域分别为211()12h h λλ++<23111()()126h h h λλλ+++<2341111()()()12624h h h h λλλλ++++<当λ为实数时，二、三、四阶显式龙格－库塔法的绝对稳定区域分别为20h λ-<<、2.510h λ-<<、 2.780h λ-<<.例8.5设有初值问题21010101(0)0xy y x x y ⎧'=-≤≤⎪+⎨⎪=⎩用四阶经典龙格－库塔公式求解时，从绝对稳定性考虑，对步长h 有何限制？解对于所给的微分方程有2100,(010)1f x x y xλ∂==-<≤≤∂+在区间[0,10]上，有201010max ||max51t x x λ<<==+由于四阶经典龙格－库塔公式的绝对稳定区间为 2.7850h λ-<<，则步长h 应满足00.557h <<.。
第8章常微分方程边值问题的数值解法
第8章常微分方程边值问题的数值解法8.1 引言第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。
只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为例介绍常用的数值方法。
一般的二阶常微分方程边值问题(boundary-value problems for second-order ordinary differential equations)为, (8.1.1)其边界条件为下列三种情况之一：(1) 第一类边界条件 (the first-type boundary conditions)：(2) 第二类边界条件 (the second-type boundary conditions)：(3) 第三类边界条件 (the third-type boundary conditions)：定理8.1.1 设(8.1.1)中的函数及其偏导数在上连续. 若(1) 对所有，有；(2) 存在常数，对所有，有,则边值问题(8.1.1)有唯一解。
推论若线性边值问题(8.1.2)满足(1)和上连续；(2) 在上,,则边值问题(8.1.1)有唯一解。
求边值问题的近似解，有三类基本方法：(1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解；(2) 有限元法(finite element method)；(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。
8.2 差分法8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法设二阶线性常微分方程的边值问题为其中在上连续，且用差分法解微分方程边值问题的过程是：(i) 把求解区间分成若干个等距或不等距的小区间，称之为单元；(ii) 构造逼近微分方程边值问题的差分格式. 构造差分格式的方法有差分法, 积分插值法及变分插值法；本节采用差分法构造差分格式；(iii) 讨论差分解存在的唯一性、收敛性及稳定性；最后求解差分方程.现在来建立相应于二阶线性常微分方程的边值问题(8.2.1), (8.2.2)的差分方程.( i ) 把区间等分，即得到区间的一个网格剖分：,其中分点，并称之为网格节点(grid nodes)；步长.( ii ) 将二阶常微分方程(8.2.2)在节点处离散化：在内部节点处用数值微分公式(8.2.3)代替方程(8.2.2)中，得, (8.2.4)其中.当充分小时，略去式(8.2.4)中的，便得到方程(8.2.1)的近似方程, (8.2.5)其中,分别是的近似值, 称式(8.2.5)为差分方程(difference equation)，而称为差分方程(8.2.5)逼近方程(8.2.2)的截断误差(truncation error). 边界条件(8.7.2)写成(8.2.6)于是方程(8.2.5), (8.2.6)合在一起就是关于个未知量，以及个方程式的线性方程组：(8.2.7)这个方程组就称为逼近边值问题(8.2.1), (8.2.2)的差分方程组(system of difference equations)或差分格式(difference scheme)，写成矩阵形式. (8.2.8)用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.7)或(8.2.8), 其解称为边值问题(8.2.1), (8.2.2)的差分解(difference solution). 由于(8.2.5)是用二阶中心差商代替方程(8.2.1)中的二阶微商得到的，所以也称式(8.2.7)为中心差分格式(centered-difference scheme).( iii ) 讨论差分方程组(8.2.7)或(8.2.8)的解是否收敛到边值问题(8.2.1), (8.2.2)的解，估计误差.对于差分方程组(8.2.7)，我们自然关心它是否有唯一解；此外，当网格无限加密，或当时，差分解是否收敛到微分方程的解. 为此介绍下列极值原理：定理8.2.1 (极值原理) 设是给定的一组不全相等的数，设. (8.2.9)(1) 若, 则中非负的最大值只能是或；(2) 若, 则中非正的最小值只能是或.证只证(1)的情形，而(2)的情形可类似证明.用反证法. 记，假设, 且在中达到. 因为不全相等，所以总可以找到某个，使，而和中至少有一个是小于的. 此时因为，所以, 这与假设矛盾，故只能是或. 证毕！推论差分方程组(8.2.7)或(8.2.8)的解存在且唯一.证明只要证明齐次方程组(8.2.10)只有零解就可以了. 由定理8.7.1知，上述齐次方程组的解的非负的最大值和非正的最小值只能是或. 而，于是证毕！利用定理8.2.1还可以证明差分解的收敛性及误差估计. 这里只给出结果：定理8.2.2 设是差分方程组(8.2.7)的解，而是边值问题(8.2.1), (8.2.2)的解在上的值，其中. 则有(8.2.11)其中.显然当时，. 这表明当时，差分方程组(8.2.7)或(8.2.8)的解收敛到原边值问题(8.7.1), (8.7.2)的解.例8.2.1 取步长，用差分法解边值问题并将结果与精确解进行比较.解因为，, 由式(8.2.7)得差分格式,, 其结果列于表8.2.1.表8.2.1准确值0 1 0 01 0.1 －0. 0332923 －0.03336562 0.2 －0. 0649163 －0.06506043 0.3 －0. 0931369 －0.09334614 0.4 －0. 1160831 －0.11634825 0.5 －0. 1316725 －0.13197966 0.6 －0. 1375288 －0.13785787 0.7 －0. 1308863 －0.13120878 0.8 －0. 1084793 －0.10875539 0.9 －0. 0664114 －0.066586510 1.0 0 0从表8.2.1可以看出, 差分方法的计算结果的精度还是比较高的. 若要得到更精确的数值解，可用缩小步长的方法来实现.8.2.2 一般二阶线性常微分方程边值问题的差分法对一般的二阶微分方程边值问题(8.2.12)假定其解存在唯一.为求解的近似值，类似于前面的做法，( i ) 把区间等分，即得到区间的一个网格剖分：,其中分点，步长.( ii ) 对式(8.2.12)中的二阶导数仍用数值微分公式代替，而对一阶导数，为了保证略去的逼近误差为，则用3点数值微分公式；另外为了保证内插，在2个端点所用的3点数值微分公式与内网格点所用的公式不同，即(8.2.13)略去误差，并用的近似值代替，，便得到差分方程组(8.2.14)其中,是的近似值. 整理得(8.2.15)解差分方程组(8.2.15)，便得边值问题(8.2.12)的差分解.特别地, 若，则式(8.2.12)中的边界条件是第一类边值条件：此时方程组(7.7.16)为(8.2.16)方程组(8.2.16)是三对角方程组，用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.16)，便得边值问题(8.2.12)的差分解.( iii ) 讨论差分方程组(8.2.16)的解是否收敛到微分方程的解，估计误差. 这里就不再详细介绍.例8.2.2 取步长，用差分法求下列边值问题的近似解，并将结果与精确解进行比较.精确解是.解因为，, 由式(8.2.17)得差分格式,, 其结果列于表8.2.2.表8.2.2准确值0 0 -0.3 -0.31 /16 -0.3137967 -0.31374462-0.3154982 -0.3154322 2/163-0.3050494 -0.3049979 3/1644-0.2828621 -0.2828427/1655-0.2497999 -0.2498180/1666-0.2071465 -0.2071930/167-0.1565577 -0.15660567/168 /2 -0.1000000 -0.10000008.3 有限元法有限元法(finite element method)是求解微分方程定解问题的有效方法之一，它特别适用在几何、物理上比较复杂的问题. 有限元法首先成功地应用于结构力学和固体力学，以后又应用于流体力学、物理学和其他工程科学. 为简明起见，本节以线性两点边值问题为例介绍有限元法.考虑线性两点边值问题其中,.此微分方程描述了长度为的可变交叉截面(表示为)的横梁在应力和下的偏差.8.3.1 等价性定理记, 引进积分. (8.3.3)任取，就有一个积分值与之对应，因此是一个泛函(functional)，即函数的函数. 因为这里是的二次函数，因此称为二次泛函.对泛函(8.3.3)有如下变分问题(variation problem)：求函数，使得对任意, 均有, (8.3.4) 即在处达到极小, 并称为变分问题(8.3.4)的解.可以证明：定理8.3.1(等价性定理)是边值问题(8.3.1), (8.3.2)的解的充分必要条件是使泛函在上达到极小，即是变分问题(8.3.4)在上的解.证 (充分性) 设是变分问题的解；即使泛函在上达到极小，证明必是边值问题(8.3.1), (8.3.2)的解.设是任意一个满足的函数，则函数，其中为参数. 因为使得达到极小，所以，即积分作为的函数，在处取极小值，故. (8.3.5)计算上式，得利用分部积分法计算积分代入式(8.3.6)，得因为是任意函数，所以必有. (8.3.8) 否则，若在上某点处有，不妨设，则由函数的连续性知，在包含的某一区间上有.作显然，且，但，这与式(8.3.7)矛盾. 于是式(8.3.8)成立，即变分问题(8.3.4)的解满足微分方程(8.3.1), 且故它是边值问题(8.3.1), (8.3.2)的解.。