第10章_常微分方程数值解

迭代法的收敛条件：
计
算 方
|
y y (k 1) i 1
(k ) i 1
|
h| 2
f
( xi1,
y(k) i 1
)
f
( xi1,
y (k 1) i 1
)
|
法 课
h 2
L|
y(k) i 1
y (k 1) i 1
|
件
其中L为Lipschitz常数。
当
0 h L 1
2
时迭代过程收敛。即迭代过程中误差没有扩大。
称为欧拉中点公式.容易看出,中点公式计算yi+1时,不仅需要yi的 值，还需要yi-1的值.
8
结束
凡是从已知(或已算出)的y0,y1,…,yi能直接从公式算出yi+1的公 式称显式，否则称隐式公式(8.4),(8.6)都是显式，(8.5)是隐式
在计算yi+1时，只需要yi的值，则称公式为单步法；若除yi之外 还需要以前的yi-1等多个值，则称多步法公式.(8.4)、(8.5)是单 计 步法，(8.6)是多步法，确切说是二步法.
0.2 xi yi
计 算
yi
1
yi
0.05
yi
2xi yi
y(0) i 1
2 xi 1 y (0)
i 1
方
y0 1
法 课
件 结果（略）
15
结束
§8.3 龙格-库塔(Runge-Kutta)法
欧拉方法是显式的一步法,使用方便,但精度较低.本节将构造出
高精度的显式一步法：龙格-库塔法,简称R-K法. 计 算 8.3.1 二阶R-K法
8.2.1 欧拉方法的导出
把区间［a,b］分为n个小区间，取步长h=(b-a)/ n ，节点 xi=x0+ih,i=0,1,2,…，n，其中x0=a，又设y (x)为上述问题的解.
5
结束
把(8.1)式两端积分，积分区间取［xi, xi+1］
xi 1
xi 1
ydx f (x, y)dx y(xi1 ) y(xi )
方 欧拉法的公式为：
法 课
yi+1=yi+h f (xi ,yi) i=0,1,2,…,n-1
件 决定其精度的是函数f (xi ,yi).如能改进这个函数,就可能提高公 式的精度.为此把公式改写成：
yi+1=yi+h (xi ,yi,h) i=0,1,2,…,n-1
(8.10)
选择函数 (xi ,yi,h),一种方法是用若干个点的函数值的线性组 合代替 (xi ,yi,h),如：
16
结束
(
xi
,
yi
,
h)
p
cjK j
j 1
K1 f (xi , yi )
计 算
j 1
K j f (xi a jh, yi h bjl Kl ) j 2,3, , p
方
l 1
法 其中cj,aj,bjl 是待定参数, aj和bjl满足
课
件
j 1
a j
b jl
j 2,3, , p
i 0,1, ,9.
y0 1
欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式
yi 1 y0 1
yi 1
2hyi
13
i 1,2, ,9.
结束
本题精确解为y=e-x,计算结果（略）
例2 用欧拉公式和梯形公式的预估校正法计算：
y y 2 x
0 x 1
计 算
y
y(0) 1
y (k 1) i 1
yi
h 2
f (xi , yi )
f
( xi1 ,
y (k) i 1
)
k 0,1,2, , i 0,1, , n 1
(8.9)
11
结束
它称为梯形公式的预估校正法，用低精度的显式欧拉公式 给出下一步yi+1的初值y(0)i+1，再用较高精度的隐式方法梯形 公式进行迭代,一般迭代一两次就能达到精度要求.
件
c1 y(xi ) c2 f (xi a2h, y(xi ) a2hy(xi )) (8.11)
由二元函数的泰勒展开：
f (xi a2h, y(xi ) a2hy(xi )) f (xi , y(xi )) a2hf x a2hy(xi ) f y O(h2 )
其中所有的偏导数都是它们在点(xi,y(xi))的值,下同 18
课 代入(8.10)
件
yi1 y(xi ) c1 c2 hy(xi ) a2c2h2 y(xi ) O(h3 )
而Taylor展开式
y( xi1 )
y( xi )
hy( xi
法能通过计算机便捷地实现，因此近年来得到迅速的发展和 广泛的应用常微分方程数值解法的特点是：对求解区间进行 剖分，然后把微分方程离散成在节点上的近似公式或近似方 程，最后结合定解条件求出近似解.因此数值解法得到的近似 解是一个离散的函数表.
1
结束
8.1.1 基本知识复习
一阶常微方程的初值问题：
y f ( x, y) a x b (1)
法 把(8.3)式积分区间改为[xi-1,xi+1]，并对右端用中矩形公式代替 课 积分，可得到：
件
y( xi1)
y(xi1) 2hf (xi , y(xi ))
h3 3
y( )
( xi , xi1)
推出：
yi1 yi1 2hf (xi , yi ) i 1,2, , n 1, (8.6)
结束
8.2.4 梯形公式及其预估校正法
对积分式(8.3)的右端用梯形积分公式近似代替，得
计
h y(xi1) y(xi ) 2
f
(xi , y(xi ))
f(xi1, y(Fra biblioteki1))h3
12
y( )
(xi , xi1)
算 方 法 课
把y(xi+1)记为yi+1，并略去余项-h3y(ξ)/12得到梯形公式：
4
结束
§8.2 欧拉(Euler)方法
对问题
计 算 方
y f (x, y)
y(a)
y0
a xb
(8.1) (8.2)
法 最简单而直观的方法是欧拉方法.欧拉方法在精度要求不高时，
课 仍不失为一实用方法.同时通过欧拉方法的讨论，容易弄清常微
件 方程初值问题数值解法的一些基本概念和构造方法的思路.
i 0,1, , n 1.
舍去Ri，并令yi=y(xi),得到：
yi1 yi hf ( xi , yi ) 这就是欧拉公式.
i 0,1, , n 1
(8.4)
6
结束
欧拉公式还可由泰勒展开得到.假定y (x)二阶连续可导，把
y(xi+1)在xi点展开：
计
y( xi1 )
y(xi )
y(xi )h
方 法 课
的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较. 方程的解析解为: y 1 2 x
件 解: 本例中欧拉公式为：
yi
1
yi
0.1
yi
2 xi yi
1.1yi
0.2 xi yi
y0 1
i 0,1, ,9.
14
结束
梯形公式只校正一次的格式为
y(0) i 1
1.1yi
a xb
y0 (1, ,m )T
3
结束
还有高阶常微分方程定解问题，如二阶定解问题：
y f ( x, y, y)
y(a
)
计 算
y(a)
方 可化作一阶常微方程组的定解问题：
法
课 件
y u( x, y) u f ( x, y, u)
y(a) , u(a)
本章着重讨论一阶常微方程初值问题的各种数值解法.这 些解法都可以写成向量形式而用于一阶常微分方程组的初值 问题.这样，也就解决了高阶方程的定解问题.
(x1, y1)，过(x1, y1)作过此点的积分曲线的切线，又与x=x2交于 (x2, y2)，…如此下去，得到一条折线，欧拉方法就是用这条折 线近似地代替曲线y(x)，故欧拉方法有时也称欧拉折线法.
8.2.2 欧拉隐式公式和欧拉中点公式
对(8.3)式的右端用右矩形公式代替积分，可得到
7
结束
y( xi1)
l 1
以上方法称为p级R-K法,选择cj,aj和bjl,可能使以上方法为p阶 方法.显然欧拉法就是一阶R-K法.
17
结束
二级R-K法的形式是:
Kyi11
yi f(
hc1K1
xi , yi )
c2
K
2
计 算
K2 f ( xi a2h, yi a2hK1)
方 此时
法 课
(xi , yi , h) c1 f (xi , y(xi )) c2 f (xi a2h, y(xi ) a2hf (xi , y(xi )))
y( xi
)
hf
( xi1,
y( xi1))
h2 2
y( )
推出：
(xi , xi1)
yi1 yi hf (xi1 , yi1 ) i 0,1, , n 1, (8.5)
计 算 方
不难看出，(8.5)式中公式两端都含有yi+1,一般情况不能由yi的值 计算yi+1,而需要解方程，故称为欧拉隐式公式.
y( )
2!
h2
算 舍去h2项，并令yi= y(xi) ，注意到y’(xi)=f (xi, yi),也得到
方 法
yi1 yi hf ( xi , yi ) i 0,1, , n 1 (8.4)
课 欧拉公式有明显的几何意义，如图8.1，过点(x0, y0)的曲线是解 件 y(x).欧拉方法是在(x0, y0)作y(x)的切线，它与直线x=x1交于
(8.3)
计 算 方 法
xi
xi
记g (x)=f (x, y (x))对(8.3) 用左矩形公式(7.3)
b g(x)dx (b a)g(a) g( ) (b a)2
a
2
课 得到： y( xi1) y( xi ) hf ( xi , y( xi )) Ri
件
Ri
y( ) h2
2
(xi , xi1)
2
结束
8.1.2 一阶常微分方程组和高阶常微分方程 一阶常微分方程组常表述为：
y1 f1 ( x, y1, , ym )
计
算 方 法 课
ym y1 (a)
fm ( x,
1
y1 ,
,
ym )
(a x b)
件
ym (a) m
写成向量形式为
y f (x, y)
y(a)
y0 ,
h
yi1 yi 2 f (xi , yi ) f (xi1 , yi1 )
(8.8)
件 它是单步隐式方法，也是一个二阶方法. 又是隐式方法，一般
求yi+1需要解方程，常用迭代法解此方程，初值y(0)i+1由显式的 欧拉公式给出,于是完整的公式写为:
y(0) i 1
yi
hf
(xi ,
yi )
算 方
8.2.3 局部截断误差和方法的阶
法 课 件
以上三个公式都有误差.假定yi,yi-1的值没有误差，用公式计算 yi+1时产生的误差称为局部截断误差，也称为单步截断误差,假 定yi+1以前的y值都没有误差，局部截断误差表示为：
R y(xi1) yi1
对欧拉公式：
R1
y( ) h2
2
O(h2 )
第8章 常微分方程的数值解
§8.1 引 言
计 科学研究和工程技术中的问题往往归结为求某个常微分方
算 程的定解问题.
方 法 课 件
常微分方程的理论指出，很多方程的定解虽然存在，但可 能十分复杂难于计算，也可能不能用简单的初等函数表示,因 此常求其能满足精度要求的近似解.
常微分方程的数值解法常用来求近似解，由于它提供的算
12
结束
例1 用欧拉方法，隐式欧拉方法和欧拉中点公式计算
y y
y(0)
1
0 x 1
计 的近似解，取步长h=0.1，并与精确值比较
算 解:欧拉方法的算式为：
方 法 课
yi 1 y0 1
yi
(1
h)
i 0,1, ,9.
件 欧拉隐式方法在本题可解出方程，不必迭代，公式为：
yi
1
yi 1 h
9
结束
对欧拉隐式公式：
R2
y( ) h2
2
O(h2 )
计 对欧拉中点公式的局部截断误差：
算
方 法 课
R3
y( ) h3
3
O(h3)
件 凡是局部截断误差为O(hp+1)的方法称为p阶方法，这样欧拉方
法和欧拉隐式方法称为一阶方法，而欧拉中点方法称为二阶
方法.一般讲高阶方法比低阶方法要精确.
10
结束
又由于： y(x) f (x, y(x))
y(x) f x f yy(x)
所以
计 f (xi a2h, y(xi ) a2hy(xi )) y(xi ) a2hy(xi ) O(h2 )
算 代入(8.11)
方 法
(xi , yi , h) c1 c2 y(xi ) a2hy(xi ) O(h2 )
计 算
y(a)
y0
(2)
方 其中f (x,y)是已知函数，(2)是定解条件，常微分方程的理论指
法 出：当f (x,y)定义在区域G=(a≤x≤b,｜y｜＜∞)，若存在正的常
课 数L，使：
件
｜f (x,y1) – f (x,y2) ｜≤L｜y1- y2|
对任意的x∈[a,b]和任意的y1, y2成立，则称f (x,y)对y满足李普 希兹(Lipschitz)条件,L称李普希兹常数. 若f (x,y)在区域G连续， 关于y满足李普希兹条件，则上述一阶常微分方程的初值问题 的解存在且唯一.我们以下的讨论，都在满足上述条件下进行.