第六章 微分方程数值解法
微分方程的数值解法
法，就产生了不同的有限差分法与不同的有限元法。
其它数学基础： 数理方程、数值代数、最优化理论与方法等
第一部分 常微分方程初值问题的数值解法
自然界与工程技术中的许多现象，其数学表 达式归结为常微分方程定解问题。 一些偏微分方程问题可以转化为常微分方程 问题来（近似）求解。 常微分方程的数值解法为偏微分方程的数值 解法提供了可供借鉴的思路。 常微分方程数值解法主要分为两大部分：
三 线性多步公式建立的基本思想
利用前面多步的信息计算 un k，以获得较 高精度的数值公式。
设 t n t0 nh ，u (t n ) 的近似值为 u n ,并记 f n f (t n , un )，k步线性多步方法一般形式为
0un 1un 1 k un k h( 0 f n 1 f n 1 k f n k )
2 ˆ ˆ ˆ ˆ ˆ1 , x ˆ2 ) 若 4 A( x1 , x2 )C ( x1 , x2 ) B ( x 为椭圆型偏微分方程
二阶偏微分方程的基本分类方法，可以推广到含 有两个以上自变量的非线性高阶偏微分方程。
初值问题的数值方法 边值问题的数值方法
这里只介绍初值问题。
目的：建立一阶常微分方程初值问题的数值解法。 模型
du f (t , u (t )) dt u (t0 ) u0 t0 t T
设初值问题的解析解 (理论解) 用 u (t n ) 表 示，数值解法的精确解用 u n 表示。其中n=1, 2, ，t n t n 1 hn 。 常微分方程初值问题的数值解是求上述初 值问题的解u(t)在点列 t n t n 1 hn 上的近似值 u n (n 0, 1, ) 。 以下设 hn不变，记为h。
计算物理学(刘金远)课后习题答案第6章：偏微分方程数值解法
第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的根是单根。
常微分方程数值解法_OK
y(xi )
O(h3)][yi
hf
(xi ,
yi )]
h2 2
y(xi ) O(h3 )
O(h2 )
欧拉法具有 1 阶精度。4
2. 隐式 Euler法
用向后差商公式代替导数项
y(xi1 ) h
y(xi )
y' (xi1 )
h 2
y' ' ( i
)
y(xi1 ) h
y(xi )
f (xi1, y(xi1 ))
i1 y(xi1 ) yi1 O(h3f)x ( x, y) f y ( x, y) f ( x, y) Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开
K2 f (xi ph, yi phK1)
f (xi , yi ) phfx (xi , yi ) phK1 f y (xi , yi ) O(h2 ) y(xi ) phy(xi ) O(h 2 )
f
(
xi
1
,
y(
xi
1
))]
h3 12
f
''( )
所以，有格式为：
yi1
yi
h[ f 2
(xi , yi )
f
(xi1, yi1 )]
上式称为梯形格式。
类似，可以算出梯形格式的误差估计式：
i1 O(h3 )
2阶的方法
梯形法是二阶、隐式单步的方法，要用迭代法求解。怎么求？
8
改进欧拉格式 /* modified Euler’s Formula */
xi1, yi h f ( xi , yi )
(i 0, ..., n 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 ，那么不但有可能提高精度，而且大大减少了计算量，这就是构造所谓线性多步法的基本思想。
第6章常微分方程数值解法
对y′=f（x，y）等式两边在（xi，xi+1）上取积分
xi1y'dx xi1 f(x,y(x)d ) x
xi
xi
《 计 算
y(xi1)y(xi)xxii1f(x,y(x)d ) x
方 法
(xi+1,yi+1)
》
(xi,yi)
xi
xi+1
第6章 常微分方程数值解法
xi 1
方 法
i0,1,2,...
》
(xi+1,yi+1)
这样得到的点列仍
(xi,yi)
为一折线，只是用 平均斜率来代替原
来一点处的斜率。
xi
xi+1
改进的欧拉公式，又称为梯形公式
第6章 常微分方程数值解法
《 计 算 方 法 》
第6章 常微分方程数值解法
不难发现:
欧拉公式 yi+1=yi +hf(x i,yi)
《 初值问题（6.1）(6.2)的数值解法的基本特点是：
计
算
求解过程顺着节点排列的顺序一步步向前推进，
方
法 》
也即按递推公式由y0，y1…..yi推yi+1，下面各种方法的
实质是建立递推公式。
第6章 常微分方程数值解法
§2 欧拉法和改进的欧拉法
一、欧拉方法
1. 基本思想
《
计 算
区间[a，b]上给定n+1个点x0，x1，x2，……xn
法
》
0.8 0.587322 0.527792
1.0 0.399383 0.367879
1.2 0.239630 0.236938
可见欧拉法的精度是很差的
第6章常微分方程初值问题的解法
ykh 2 k[ (ykx k 1 ) ( yk 1x k 1 1 )]
yk11 29 1yk1k05110
预估－校正Euler方法：
y k 1 0 .90 y k 5 0 .00 k 9 0 .1 5
20
Euler方法
xk
yk
yk y(xk)
0.0 1.000000
0.0
梯形方法
yk
yk y(xk)
1.000000
0.0
续
预估－校正方法
yk
yk y(xk)
1.000000
0.0
0.1 1.000000 0.2 1.010000
4.8×10－3 8.7×10－3
1.004762 1.018594
y(0) 1
其解析解为： y1xe-t2dt x[0,1] 0 很难得到其解析解
4
例如：
y=x+y , x[0,1]
y(0) 1
其解析解为 yx12ex
只有一些特殊类型的微分方程问题能够得到用解析表达式 表示的函数解，而大量的微分方程问题很难得到其解析解。
因此，只能依赖于数值方法去获得微分方程的数值解。
例如：
y=x+y , x[0,1]
y(0) 1
其解析解为：yx12ex
3
但是, 只有一些特殊类型的微分方程问题能够得到用解析 表达式表示的函数解，而大量的微分方程问题很难得到其解 析解。
因此，只能依赖于数值方法去获得微分方程的数值解。
例如：
y =e-x2 ,
x[0,1]
7.5×10－5 1.4×10－4
《微分方程的数值解》课件
谱方法：将微分方程离散化为谱方程， 然后求解
边界元法：将微分方程离散化为边界 元方程，然后求解
有限元法：将微分方程离散化为有限 元方程，然后求解
网格法：将微分方程离散化为网格方 程，然后求解
数值解法的步骤
确定微分方程的初值 和边界条件
选择合适的数值解法， 如欧拉法、龙格-库塔 法等
实解
应用：广泛应 用于工程、物 理、化学等领
域
优缺点：优点 是计算速度快， 缺点是精度较
低
非线性方程的数值解法
牛顿法：通过迭 代求解非线性方 程
拟牛顿法：通过 迭代求解非线性 方程，比牛顿法 收敛更快
割线法：通过迭代 求解非线性方程， 适用于求解单变量 非线性方程
迭代法：通过迭 代求解非线性方 程，适用于求解 多维非线性方程
05 数值解法的实现
M AT L A B 编 程 实 现
MATLAB简介： MATLAB是一种高 级编程语言，广泛 应用于科学计算、 数据分析等领域
数值解法：包括欧 拉法、龙格-库塔 法、四阶龙格-库 塔法等
MATLAB实现：使 用MATLAB编写程 序，实现数值解法 的计算
示例代码：给出 MATLAB实现数值 解法的示例代码， 并解释其含义和作 用
设定时间步长和空间 步长
计算微分方程的解， 并进行误差分析
绘制解的图形，并进 行结果分析
对比不同数值解法的 优缺点，选择最优解 法
04 常用的数值解法
欧拉方法
基本思想：将微分 方程转化为差分方 程，然后求解差分 方程
优点：简单易行， 适用于初值问题
缺点：精度较低， 稳定性较差
改进方法：改进欧 拉方法，如改进欧 拉方法、龙格-库 塔方法等
第六章—常微分方程的数值解法 PPT
§6.1 引言
初值问题的数 点值 ：解 按法 节特 点顺 进序 ，依 由
知的 0,yy1,,y， i 求i出 1 ，y这可以通过 得递 到推 。
初值 问题 的 常见 解法
单步法： 利用前一个单步的信息(一个点)，在y=f(x)
上找下一点yi， 有欧拉法，龙格－库格法。
预测校正法： 多步法，利用一个以上的前点信息求f(x)
第六章 常微分方程的数值解法
本章内容
§6.1 引言 §6.2 欧拉方法 §6.3 龙格—库塔方法 §6.4 边值问题的数值方法
§6.1 引言
一. 问题提出
有一个或多个导数及其函数的方程式称为微分 方程，在工程中常遇到求解微分方程的问题。
如，一阶常微分初方值程问的题 dy f(x,y) x[a,b] dx y(x0)y0
推进Pn1(xn1, yn1， ) 显然两个顶点P， n Pn1的坐标有关系
yn1 - yn xn1 - xn
f (xn, yn),
即yn1 ynr)公式。
y
y y(x)
P2 P3 P4 Pn
P1
P0
x O
§6.2 欧拉方法及其改进
例：利用 Euler 方法求初值问题
y(x0)y(x1) hy(x0)
记为
y ( x 1 ) y ( x 0 ) h y ( x 0 ) y 0 h f ( x 0 ,y 0 ) y 1
y i 1 y i h f ( x i,y i)( i 0 ,.,. n . 1 )
几何意义：折线逼近解y y(x)曲线。
设已做出折线的顶点P， n 过Pn(xn, yn)依方向场的方向再
需要用数值方法来求解，一般只要求得到若干个 点上的近似值或者解的简单的近似表达式(精度要求 满足即可)。
第六章_常微分方程初值问题的数值解法_习题课
h2 h3 y ( x n ) y ( x n ) O(h 4 ) 2 6 而且 y ( x n ) f ( x n , y ( x n )) ， y ( x n 1 ) f ( x n 1 , y ( x n 1 )) ，对 y ( x n 1 ) 也在 x n 处作 Talor 展开， y ( x n 1 ) y ( x n ) hy ( x n )
湖北民族学院理学院《数值计算方法》教学辅导材料
陈以平编写
h2 h3 y ( x n ) y ( x n ) O(h 4 ) 2 6 h h h2 h3 y ( x n ) y ( x n ) y ( x n ) y ( x n ) y ( x n ) O(h 4 ) 2 2 2 12 h3 y ( x n ) O(h 4 ) O(h 3 ) 12 h3 所以，梯形公式是 2 阶方法，其截断误差的主项是 y ( x n ) 。 12 y ( x n ) hy ( x n )
y k (0.9 0.1y k sin x k ) 0.1( y k 1 y k 1 sin x k 1 )
2
当 k=0，x0=1, y0=1 时，x1=1.2，有 y y (. . y sin x ) (. sin ) .
y f ( x, y ) 3.求解初值问题 欧拉法的局部截断误差是( y ( x ) y 改进欧拉法的局部截断误差是( ); 四阶龙格－库塔法的局部截断误差是( ). (A)O(h2) (B)O(h3) (C)O(h4) (D)O(h5)
4. 改进欧拉法的平均形式公式是( ) y p y k hf ( x k , y k ) y p y k hf ( x k , y k ) (B) y c y k hf ( x k , y p ) .(A) y c y k hf ( x k , y p ) y k ( y p y c ) y k ( y p y c ) y p y k hf ( x k , y k ) y p y k hf ( x k , y k ) (C) y c y k hf ( x k , y p ) (D) y c y k hf ( x k , y p ) y k h ( y p y c ) y k ( y p y c ) (D) 答案：
微分方程数值解法-课件
只要f ( x, y)在a,bR1上连续, 且关于 y 满足
Lipschitz 条件，即存在与 x , y 无关的常数 L 使
|f ( x ,y 1 ) f ( x ,y 2 ) | L |y 1 y 2 |
对任意定义在 a , b 上的 y1x,y2x都成பைடு நூலகம்，
（2）用代数的方法求出解函数 y解y(在的x)唯点一的x性k近似值 y k y (x k ) 工k 程 师1 ,关2 , 注,n 解yk的光滑性 解*的振动y( 性xk )
解的周期性
数学界关注
yy解(x的) 稳定性 解的混沌性
……
所谓数值解法：
求函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值
称 f ( x, y) 在区域D上对 y 满足Lipschitz条件是指:
L0 s.t. f(x,y1)f(x,y2) Ly1y2 ,
x[a,b],y1,y2[y(x),y(x)]
4、 迭代格式的构造
(1) 构造思想：将连续的微分方程及初值条件离散为线性方程 组加以求解。由于离散化的出发点不同，产生出各种不同的数 值方法。基本方法有：有限差分法（数值微分）、有限体积法 （数值积分）、有限元法（函数插值）等等。
离变量法”等特殊方法求得初等函数形式的解，绝大
部分方程至今无法理论求解。
如
y ' sx i2 n )y , y '( 1 x sy i,ny ' e x 2 xy 等等
2、数值解的思想
如果找不到解函数
（1）将连续变量 x[a离,b散] 为 数学界还关注：
a x 0 x 1 x k x n 解 的b 存在性
计算物理学：第六章 偏微分方程的数值解法
常数: a
format long; h = (maxx-minx)/(n-1); if a>0
精度： O(Δt, Δx)
差分方程的稳定性和收敛性：
收敛性：理论上，h → 0 时，解逼近准确解
稳定性：初值有小干扰的情况下，干扰不会被扩大 传播，而是被“磨灭”
对流方程
迎风格式:
∂u + a ∂u = 0 ∂t ∂x
二层显式格式
un+1 k
=
ukn
−
aΔt Δx
(uk +1
−
uk )
un+1 k
区间数： n = 1 = 50 0.02
初始值 : u0 , u1, u2 L, u50
t = 0 时，
⎧10 x + 1 − 0.1 ≤ x ≤ 0
U ( x)
=
⎪ ⎨‐
10
x
+
1
0 ≤ x ≤ 0.1
⎪⎩0
其余
⎧u0 = 1
⎪ ⎪⎪ ⎨
u1 u2
= =
−10 × 0.02 + 1 −10 × 2 × 0.02
h2
2.微分方程离散化： 差分公式:
dui = ui+1 − ui + O(h)
dx
h
dui = ui − ui−1 + O(h)
dx
h
dui = ui+1 − ui−1 + O(h2 )
dx
2h
dui = − ui+2 + 4ui+1 − 3ui + O(h2 )
dx
2h
d 2 ui dx 2
第六章常微分方程数值解法
三阶龙格-库塔法和四阶龙格-库塔法的计算公式也是有无 穷多个。
四阶龙格-库塔法计算举例 四阶龙格 库塔法计算举例
例 取步长h=0.2,用四阶龙格-库塔公式求下面初值问题的数值解。
解
y' = y − 2yx (0 ≤ x ≤1 ) y(0) =1
由标准的四阶龙格-库塔公式 0, y =1 h = 0.2 f (x, y) = y − 2x / y, x0 = 0 ,
第三节 线性多步法
前面讨论的各种方法,计算yi+1时只使用yi的值,统称为
单步法。在计算yi+1时我们可以使用前面的k个yi的值，即
由yi ， yi-1 ，yi -2，…, yi-k+1计算yi+1 ，这便称为多步法。 下面讨论常用的四阶阿当姆斯（亚当斯）线性多步法。
线性多步法的基本思想
我们知道初值问题 y' = f (x, y) y(x0 ) = y0 与积分方程
h 2
龙格-库塔法 第二节 龙格 库塔法
欧拉公式 可写为
y0 = y(x0 ) (i = 0,1 2,...) , yi+1 = yi + hf (xi , yi )
yi+1 = yi + k1 其 部 断 差 o(h2 )。 局 截 误 为 k1 = hf (xi , yi )
p yi(+1) = yi + hf (xi , yi )(预 ） 测 (c) h p yi+1 = yi + [ f (xi , yi ) + f (xi+1, yi(+1) )] 校 ） （ 正 yi + k1 + k2 2 k1 = hf (xi , yi ) 其局部 截断误 差为 (h3 ) o k = hf (x + h, y + k ) i i 1 2
第六章常微分方程的数值解法
第六章常微分方程的数值解法第六章常微分方程的数值解法在自然科学研究和工程技术领域中，常常会遇到常微分方程的求解问题。
传统的数学分析方法仅能给出一些简单的、常系数的、经典的线性方程的解析表达式，不能处理复杂的、变系数的、非线性方程，对于这些方面的问题，只能求诸于近似解法和数值解法。
而且在许多实际问题中，确确实实并不总是需要精确的解析解，往往只需获得近似的解或者解在若干个点上的数值即可。
在高等数学课程中介绍过的级数解法和逐步逼近法，能够给出解的近似表达式，这一类方法称为近似解法。
还有一类方法是通过计算机来求解微分方程的数值解，给出解在一些离散点上的近似值，这一类方法称作为数值方法。
本章主要介绍常微分方程初值问题的数值解法，包括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 -≈。
《微分方程数值解法》PPT课件
方程的解 U~n 。为了弄清差分格式（2.58）的稳定性条件， 给出稳定的定义：
对于任意给定的 0 ，存在与h, k 无关且依赖于 的
正数 ，使当
U~0 U 0 V 0
时，对于任何的 n0 nk T ，差分格式得到的解U~ n ,U n
满足不等式
U~n U n V n
连同初值条件：U
0 m
mk
, m
1,2,M
1
边值条件：U
n 0
U
n M
0, n
0,1,2,, N
逐层解出结点处的U 值。
现在对
h
ห้องสมุดไป่ตู้
，取二种 20
k
，使
r
k h2
5 和5 11 9
。图2.9
和图2.10中的曲线表示不同的时刻微分方程的精确解，图
中“ ”表示差分方程的解
（2.54）
下面我们先研究上式右边第二项，即差分方程的理论
解与计算机上解得的近似解之间的差别是随着n 的增大而
无限增大还是有所控制。如果这种差别是无限增加，则称
差分格式不稳定，显然不稳定的格式是不能使用的，因为
误差的无限增加淹没了真解。上例中的r 5 时就是差分
9
方程不稳定的情况。从差分方程，比如格式（2.29）可知，
如果差分方程为显式，则对所有的n ，An I ；如 An I
果 An
0
，
U U
n1 CnU
0
n
An1en , Cn
An1Bn
（2.58）
，则隐式格式可以写成显式形式。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
得堡科学院从事研究工作，在分析学、
数论、力学方面均有出色成就，并应 俄国政府要求，解决了不少地图学、 造船业等实际问题。 24岁晋升物理学教授。
1735年（28岁）右眼失明。
1741年 - 1766（34岁-59岁）任德国科学院物理数学所所
长，任职25年。在行星运动、刚体运动、热力学、弹道学、人 口学、微分方程、曲面微分几何等研究领域均有开创性的工作。 1766年应沙皇礼聘重回彼得堡，在1771年（64岁）左眼失 明。 Euler是数学史上最多产的数学家，平均以每年800页的速 度写出创造性论文。他去世后，人们用35年整理出他的研究成
5、微分方程的数值解法需要解决的主要问题
(1) 如何将微分方程离散化，并建立求其数值解的迭代公式？
(2) 如何估计迭代公式的局部截断误差与整体误差？ (3) 如何保证迭代公式的稳定性与收敛性?
二、初值问题解的存在唯一性
考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:
隐式欧拉法的局部截断误差：
Ri y( xi 1 ) yi 1
h y( xi ) O(h3 ) 2
2
即隐式欧拉公式具有 1 阶精度。
梯形公式 /*trapezoid formula */
— 显、隐式两种算法的平均
h yi 1 yi [ f ( xi , yi ) f ( xi 1 , yi 1 )] ( i 0, ... , n 1) 2
可用来检验近似解的准确程度。 从上表最后一列，我们看到取步长 h 0.1
进行计算，数值解已达到了一定的精度。
3、 欧拉公式的改进： 隐式欧拉法 /* implicit Euler method */
向后差商近似导数
y( x1 ) y( x1 ) y( x0 ) h
y ( x1 ) y 0 h f ( x1 , y ( x1 ))
yk y( xk ) k 1,2,, n 工程师关注
解的光滑性 yk y( xk ) 解的振动性 * 解的周期性
数学界关注
y y解的稳定性 ( x)
解的混沌性 ……
所谓数值解法：
求函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值
yi y( xi ) (i 1, ... , n)
3 R y x y O h 注：梯形公式的局部截断误差 i i 1 i 1 ，
即梯形公式具有2 阶精度，比欧拉方法有了进步。 但注意到该公式是隐式公式，计算时不得不用到
迭代法，其迭代收敛性与欧拉公式相似。
中点欧拉公式 /* midpoint formula */
中心差商近似导数
算法有p 阶精度。
欧拉法的局部截断误差：
yi y( xi )
y '( xi ) f ( xi , y( xi ))
Ri y( xቤተ መጻሕፍቲ ባይዱ 1 ) yi 1
[ y( xi ) hy( xi )
h2 2
y( xi ) O(h3 )] [ yi hf ( xi , yi )]
果74卷。
2、欧拉法的局部截断误差
定义2.1 在假设 yi = y(xi)，即第 i 步计算是精确的前提
下，考虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断 误差 /* local truncation error */。
定义2.2 若某算法的局部截断误差为O(hp+1)，则称该
dy f ( x, y) dx y ( a ) y0 x [a , b ]
只要 f ( x, y) 在a, b R1上连续, 且关于 y 满足 Lipschitz 条件， 即存在与 x, y 无关的常数 L 使
| f ( x, y1 ) f ( x, y2 ) | L | y1 y2 |
L 0
s.t.
f ( x, y1 ) f ( x, y2 ) L y1 y2 , x [a, b], y1 , y2 [ y ( x) , y ( x) ]
4、 迭代格式的构造
(1) 构造思想：将连续的微分方程及初值条件离散为线性方程 组加以求解。由于离散化的出发点不同，产生出各种不同的数 值方法。基本方法有：有限差分法（数值微分）、有限体积法 （数值积分）、有限元法（函数插值）等等。 (2) 一般构造方法： 离散点函数值集合 + 线性组合结构 → 近似公式
·····
pk
pn pn1
p0
Euler格式
yn1 yn f ( xn , yn ) xn1 xn
yn1 yn hf ( xn , yn )
x0 x1 xk xn xn 1
18世纪最杰出的数学家之一，13岁 时入读巴塞尔大学，15岁大学毕业， 16岁获得硕士学位。 1727年-1741年（20岁-34岁）在彼
第6章
§6.1 §6.2
常微分方程数值解法
引 言 欧拉方法
§6.3
龙格—库塔方法
§6.1 引 言
微分方程数值解一般可分为：常微分方程数值解和偏微分 方程数值解。自然界与工程技术中的许多现象，其数学表达式 可归结为常微分方程（组）的定解问题。一些偏微分方程问题 也可以转化为常微分方程问题来（近似）求解。Newton最早采 用数学方法研究二体问题，其中需要求解的运动方程就是常微 分方程。许多著名的数学家，如 Bernoulli（家族），Euler、 Gauss、Lagrange和Laplace等，都遵循历史传统，研究重要 的力学问题的数学模型，在这些问题中，许多是常微分方程的 求解。作为科学史上的一段佳话，海王星的发现就是通过对常 微分方程的近似计算得到的。本章主要介绍常微分方程数值解 的若干方法。
的方法称为微分方程的数值解法。
y1 ,, yn 称为微分方程的数值解。
称节点间距 hi xi 1 xi (i 0, ... , n 1) 为步长，
通常采用等距节点，即取 hi = h (常数)。
3、相关定义 记 D {( x, y) a x b, y( x) y y( x) } 称 f ( x, y) 在区域D上对 y 满足Lipschitz条件是指:
方 法 显式欧拉 隐式欧拉 梯形公式 中点公式
简单 稳定性最好
精度提高 精度提高, 显式
精度低 精度低, 计算量大
计算量大 多一个初值, 可能影响精度
改进欧拉法 /* modified Euler’s method */
Step 1: 先用显式欧拉公式作预测，算出
yn1 yn hf xn , yn
sin y
,
y' e
x 2 xy
等等
2、数值解的思想
如果找不到解函数
（1）将连续变量 x [a , b] 离散为 数学界还关注：
a x0 x1 xk xn解的存在性 b
y( x ) 在 xk 点的近似值 （2）用代数的方法求出解函数 y 解的唯一性
一、初值问题的数值解法
1、常微分方程与解
如果函数 y y( x ) 在区间[a,b]内n阶可导，称方程
F ( x, y, y' , y'' ,, y( n) ) 0
为n阶常微分方程。
满足方程的函数 y y( x ) 称为微分方程的解。 如
y' 2 x
2 y x C , (C为任意常数） 则
y' f ( x , y ) y ( x 0 ) y0
定解条件（初始条件）
是否能够找到定解问题的解取决于 f ( x , y )
仅有极少数的方程可以通过“常数变易法”、“可分
离变量法”等特殊方法求得初等函数形式的解，绝大
部分方程至今无法理论求解。
如
y' sin(xy ),
2
y' 1 x
§6.2 Euler方法 1、Euler格式
第一步：连续变量离散化
y' f ( x , y ) y ( x 0 ) y0
x x0 , x1 ,, xk , xn ,
p1
第二步：用直线步进
y1 y0 f ( x0 , y0 ) x1 x0
y2 y1 f ( x1 , y1 ) x2 x1
h2 3 y ( xi ) O(h ) 2
O(h )
2
Ri 的主项
/* leading term */
欧拉法具有 1 阶精 度。
例1: 用欧拉公式求解初值问题
2 y 2 xy y 0 1
（0 x 1.2）
取步长 h 0.1。 解: 应用Euler公式于题给初值问题的具体形式为：
y xi yi
0.009901 0.018462 0.024153 0.026320 0.025250 0.021852 0.017213 0.012262 0.007626 0.003642 0.000422 0.002053
2 这个初值问题的准确解为 y x 1 1 x ，
yi
1.000000 0.980000 0.941584 0.888389 0.825250 0.757147 0.688354 0.622018 0.560113 0.503642 0.452911 0.407783
y xi
0.990099 0.961538 0.917431 0.863069 0.800000 0.735294 0.671141 0.609756 0.552486 0.500000 0.452489 0.409836