初值问题的数值解法
常微分方程初值问题的数值解法
1 1 2 1 , 2 p 2
这里有 3 个未知 数， 2 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到，p 1, 1 2 1 就是改进的欧拉法。
Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开
K 2 f ( xi ph, yi phK1 ) f ( xi , yi ) phf x ( xi , yi ) phK1 f y ( xi , yi ) O( h2 )
y( xi ) phy( xi ) O(h2 )
d f ( x, y) dx 首先希望能确定系数 1、2、p，使得到的算法格式有 2阶 dy 精度，即在 yi y( xi ) 的前提假设下，使得 f x ( x, y) f y ( x, y) dx Ri y( xi 1 ) yi 1 O(h3 ) f x ( x, y) f y ( x, y) f ( x, y) y( x )
y( x0 ) y0 yk 1 yk h f ( xk , yk 1 )
, k 0,1,...
隐式欧拉法的求解: 利用迭代的思路进行.
yi 1 yi hf ( xi , yi 1 )
变换为
y
( k 1) i 1
yi hf ( xi , y )
y i 1 K1 K2
1 1 y i h K 1 K 2 2 2 f ( xi , yi ) f ( xi h, yi hK 1 )
步长一定是一个h 吗？
§2 Runge-Kutta Method
第7章 常微分方程初值问题的数值解法
例1 函数f ( t , y ) = t y 在区域D0 = {( t , y ) | 1 ≤ t ≤ 2, −3 ≤ y ≤ 4}
关于y满足Lipschitz条件，相应的Lipschitz常数可取为L = 2
3 存在性定理 定理1 设函数f ( t , y )在凸集D ⊂ R 2中有定义，若存在常数
(7.2.7)
称为显式Runge-Kutta(龙格-库塔 )方法，简称R-K方法，
其中正整数N 称为R-K方法的级，所有ci , ai , bij 都是待定 常数。
根据定义(7.2.7)，N 级R-K方法(7.9)的局部截断误差为
Rn+1 = y( t n+1 ) − y( t n ) − h∑ ci ki
dy 其斜率为 = f ( t0 , y0 ) dt ( t0 , y0 ) 由 点 斜 式 写 出 切线 方 程 dy y = y0 + ( t − t0 ) = y0 + ( t − t0 ) f ( t0 , y0 ) dt ( x0 , y )
0
等步长为h,则t1 - t0 = h, 可由切线算出 y1 ： 则 y1 = y0 + hf ( t0 , y0 ) 按此逐步计算y( tn ), 在tn +1处的值 ： yn+1 = yn + hf ( tn , yn ) y 注意： 这是“ 注意 ： 这是 “ 折 yN 线法” 而非“ 线法 ” 而非 “ 切 线法” 线法 ” 除第一个 点是曲线切线外， 点是曲线切线外 ， 其他点不是切线 y2 而是折线(如右 y1 y0 图所示)。 图所示 。
பைடு நூலகம்
则称数值解法(7.5)为显式方法。否则，称数值解法(7.3) 为隐式方法。
第九章_常微分方程初值问题的数值解法
向与已给向量场在该点的方向一致。
而欧拉递推公式 (9.7) 所定义的点 ( xi , yi ) 可以看作 是通过准确初始点 ( x0 , y0 ) 的曲线多边形的顶点，而 每根连线的方向由它左边终点的向量场决定。
y
y y ( x)
o
x0
x1
x2
图9.2
xn
x
由于欧拉方法是用一条折线近似地代替解 y ( x) ， 所以欧拉方法也叫做折线法。
(9.23) , n 1
称为 r 阶泰勒( Taylor )方法。当 r 1 时就是欧拉方法的 公式，现取 r 2 ，则有
y0 y( x0 ) 2 h yi1 yi hf ( xi , yi ) ( f x( xi , yi ) f ( xi , yi ) f y( xi , yi )) i 0,1,2, , n 1 2
令
i1 yi1 yi*1
(9.17)
* i 1
* * y y 其中, i 1 是 yi 1 的近似值，若 i 1 不大，则 i 1 可以
作为 y( xi1 ) 的近似值， 若 i 1 很大，y 就不能用。
定义用一个数值方法，求解微分方程
y y
其中, 是一个复常数， 对给定步长 h 0 ， 在计算 yi 时 引入了误差 i 。 若这个误差在计算后面的 yi k (k 1, 2, ) 中所带来的误差按绝对值均不增加， 就说这个数值方 法对于这个步长 h 和复数 是绝对稳定的。
(9.26)
(9.26)式中的 yi 1 可以用欧拉方法得到，即得
(0) y i 1 yi hf ( xi , yi ) ( k 1) (k ) y y hf ( x , y i i 1 i 1 ) i 1
数值分析课件第9章常微分方程初值问题数值解法
上页 下页
设用欧拉公式
y(0) n1
yn
hf
( xn ,
yn )
给出迭代初值 yn(0)1，用它代入(2.5)式的右端，使之转
化为显式，直接计算得
y(1) n1
yn
hf
( xn1 ,
y(0) n1
),
然后再用 yn(代1)1 入(2.5)式，又有
在yn=y(xn)的前提下，f(xn,yn )=f(xn,y(xn))=y(xn).于是
可得欧拉法(2.1)的公式误差为
y(xn1)
yn1
h2 2
y(n )
h2 2
y(xn ),
(2.3)
称为此方法的局部截断误差.
上页 下页
如果对方程(1.1)从xn到xn+1积分，得
y( xn1) y( xn )
基于上述几何解释，我们从初始点P0(x0, y0)出发, 先依方向场在该点的方向推进到x=x1上一点P1，然后 再从P1点依方向场在该点的方向推进到 x=x2 上一点 P2 , 循环前进做出一条折线P0 P1 P2.
上页 下页
一般地，设已做出该折线的顶点Pn,过Pn(xn, yn)依
方向场的方向再推进到Pn+1(xn+1, yn+1)，显然两个顶
(2.13)
校正 yn1 yn 2 [ f ( xn , yn ) f ( xn1, yn1 )] .
或表为下列平均化形式
yp yc
y(2) n1
yn
hf
( xn1 ,
y(1) n1
).
如此反复进行，得
y(k1) n1
微分方程初值问题的数值解法
积分法:
yk 1 yk h f ( xk , yk ) y ( x0 ) y0
积分项利用矩形公式计算
(1) y( xk 1 ) y( xk )
xk 1
xk
f (t , y(t ))dt
(★)
xk 1
xk
f (t , y(t ))dt h f ( xk , yk ) y( xk 1 ) y( xk ) h f ( xk , yk )
引言
初值问题的数值解法:求初值问题的解在一系列节点的值 y ( xn )的近似值 yn 的方法.本章数值解法的特点:都是采用“步进 式”,即求解过程顺着节点排列的次序一步步向前推进. 常微分方程初值问题: dy f ( x, y ), x [a, b] dx y ( x0 ) y0
替 f (x , y)关于 y 满足Lipschitz条件. 除了要保证(1)有唯一解外,还需保证微分方程本身是稳定的,即 (1)的解连续依赖于初始值和函数 f (x , y). 也就是说, 当初始值 y0 及函数 f (x , y)有微小变化时, 只能引起解的微小变化.
注: 如无特别说明,总假设(1)的解存在唯一且足够光滑. 在 f 连续有界, 则 f (x , y)对变量 y 可微的情形下, 若偏导数 y 可取L为
也称折线法 x
2. 梯形法
若采用梯形公式计算(★)中的积分项,则有 h y ( xk 1 ) y ( xk ) [ f ( xk , y ( xk )) f ( xk 1 , y ( xk 1 ))] 2 h yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk 1 )] 2 称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做 法: (0) (0) 先用Euler法求出初值 yk ,1 即 ,将其代入梯形公式 yk 1 yk h f ( xk , yk ) 的右端,使之转化为显式公式,即 h ( l 1) (l ) yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk (☆ ) 1 )] 2
第九章常微分方程初值问题的数值解法
当r=0时,称为单步法,这恰好就是欧拉方法.
3 当r =1时,称为二步法,其局部截断误差为 O h .
Adams线性多步法的特点:
r r , 1. 对任意的 j 1
r j 0
2. Adams线性多步法的截断误差的来源是数值积分的截断误差, 而数值积分的截断误差来源于插值误差,因此阿达姆斯方法的 的局部截断误差是 O hr 2 3. Adams线性多步法与Runge-Kutta法比较,R-K每提高一阶精 度至少需要多计算一个函数值;而线性多步法每提高一阶精 度只需多用一个已知的数据,所以从这个角度看计算量较小, 但线性多步法必须在单步法的基础上做. 所以在实际计算时 常常把两者相结合. 13
i
再对右端的积分应用数值积分(用函数值的线性组合来近似积分)
r yi 1 yi h j f xi j , yi j r j 0
12
由于 j 仅与r , j 有关, 而与f , i 无关，所以可以事先制成表。
r
当r =2时,称为三步法,其局部截断误差为 O h
i i
这个方法称二阶泰勒方法
6
2. 梯形方法和改进欧拉方法
xi 1 x i
h y x dx y xi 1 y xi [ f xi , yi f xi 1, yi 1 ] 2 h yi 1 yi [ f xi , yi f xi 1 , yi 1 ] 2
常微分方程初值问题解的存在唯一性定理:
以后我们总假定给出的方程都满足该定理的条件. 微分方程初值问题数值解法的特点： 先把方程离散化,即在区间[a,b ]中插入一些节点(通常采用等 距节点) ba
第9章 常微分方程初值问题数值解法
oa
b
a f ( x)dx (b a) f (b)
中矩形公式
b
ab
a f ( x)dx (b a) f ( 2 )
计算方法
梯形公式
bx
右矩形公式 中矩形公式 左矩形公式
§ 欧拉方法几何意义
y y y(x)
y0 y1 y2 0 x0 x1 x2
计算方法
x
§ 隐式欧拉方法
➢隐式欧拉法 /* implicit Euler method */
初 值 问 题 的 解 必 存 在 且唯 一 。
计算方法
§9.1 引言
三. 数值解法含义
所谓数值解法, 就是设法将常微分方程离散化, 建 立差分方程, 给出解在一些离散点上的近似值。
微分方程的数值解: 设方程问题的解y(x)的存在区 间是[a,b], 令a= x0< x1<…< xn =b, 其中hk=xk+1-xk, 如是等距节点h=(b-a)/n, h称为步长。
yi1 yi1 2h f ( xi , yi ) i 1, ... , n 1
计算方法
预估-校正法
三. 预估 — 校正法
/* predictor-corrector method */
方法 显式欧拉 隐式欧拉 梯形公式
中点公式
简单
稳定性最好
精度提高
精度低
精度低, 计算量大
计算量大
精度提高, 显式
在x0 x X上的数值解法。
四. 误差估计、收敛性
和稳定性
计算方法
§9.2 简单的数值方法与基本概念
一. 欧拉（Euler）格式
设 节 点 为xi a ih (i 0,1,2 , n) 方 法 一 ：Taylor展 开 法
专业计算方法第九章-常微分方程初值问题的数值解法-2012.2
第九章常微分方程初值问题的数值解法§1 相关知识§2 几种简单的数值方法§3 Runge–Kutta方法§4 单步法的进一步讨论§5 线性多步法西北工业大学理学院欧阳洁1上述定理称为一阶常微分方程初值问题解的适定性（存在性、惟一性与稳定性）定理。
对所讨论的一阶常微分方程初值问题，本章假设该问题是适定的，即解析解y(x)在区间[a,b]上是存在、惟一，且具有充分的光滑度。
因此f(x,y(x))也充分光滑。
西北工业大学理学院欧阳洁3西北工业大学理学院欧阳洁5常微分方程初值问题的数值解法分为：①单(一)步法：计算时，只用到和，即前一步的值。
1+n y n y n n x x ,1+显式单步法的一般形式为②多步法：计算时，除用到和以外，还用到和，即用到前k 步的值。
p n x −)1;1,2,1(>−=−k k p y p n L 1+n y n y n n x x ,1+对单步法与多步法，有显式与隐式方法之分。
显式、隐式多步法的一般形式类似。
隐式单步法的一般形式为),,(1h y x h y y n n n n ϕ+=+),,,,(111h y x y x h y y n n n n n n ++++=ϕ数值解法建立的过程：通过一定的离散化方法,将连续性问题的求解转化为有限个离散节点上解析解近似值的求解。
常用的离散化方法：Taylor 展开法；差商直接代替微商；数值积分法。
西北工业大学理学院欧阳洁16设一般的单步法为：显式公式隐式公式定义为某一数值方法在处的整体截断误差。
111()n n n e y x y +++=−1n x +五单步法的局部截断误差和阶整体截断误差不仅与这一步的计算有关，还依赖于前面各步的计算。
1+n x 011,,,,x x x x n n L −下面着重分析计算中某一步的误差（局部截断误差），后面将对显式单步法，给出整体截断误差与局部截断误差之间的关系。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第六章 常微分方程初值问题的数值解法_习题课1 .欧拉法的局部截断误差的阶为 。
改进欧拉法的局部截断误差的阶为 。
三阶龙格-库塔法的局部截断误差的阶为 。
四阶龙格-库塔法的局部截断误差的阶为 。
2. 欧拉法的绝对稳定实区域为 。
二阶龙格-库塔法的绝对稳定实区域为 。
三阶龙格-库塔法的绝对稳定实区域为 。
四阶龙格-库塔法的绝对稳定实区域为 。
3.求解初值问题欧拉法的局部截断误差是( );⎩⎨⎧=='00y x y y x f y )(),(改进欧拉法的局部截断误差是( );四阶龙格－库塔法的局部截断误差是( ). (A)O (h 2) (B)O (h 3) (C)O (h 4) (D)O (h 5)4. 改进欧拉法的平均形式公式是( ).(A)⎪⎪⎪⎩⎪⎪⎪⎨⎧+21=+=+=1+)(),(),(c p k p k k c k k k p y y y y x hf y y y x hf y y (B)⎪⎪⎪⎩⎪⎪⎪⎨⎧+21=+=+=1+1+1+)(),(),(c p k p k k c k k k p y y y y x hf y y y x hf y y(C)⎪⎪⎪⎩⎪⎪⎪⎨⎧+2=+=+=1+1+)(),(),(c p k p k k c k k k p y y h y y x hf y y y x hf y y (D)⎪⎪⎪⎩⎪⎪⎪⎨⎧+21=+=+=1+1+)(),(),(c p k p k k c k k k p y y y y x hf y y y x hf y y 答案：(D)5. 解微分方程初值问题的方法，( )的局部截断误差为O (h 3). (A) 欧拉法 (B)改进欧拉法 (C)三阶龙格－库塔法 (D) 四阶龙格－库塔法 答案：(B)解答：改进欧拉法的局部截断误差是二阶精度，O(h3)。
6. 对Euler 公式推导局部截断误差及其主项,并指出该方法是几阶方法。
解：其局部截断为))(,()()(11n n n n n x y x hf x y x y T --=++ 对在处作Talor 展开，有)(1+n x y n x )()(2)()()(321h O x y h x y h x y x y n n n n +''+'+=+而且，因此其局部截断为 ))(,()(n n n x y x f x y =' ))(,()()(11n n n n n x y x hf x y x y T --=++)()(2)()(32h O x y h x y h x y n n n +''+'+=)()(n n x y h x y '--)()(232h O x y h n +''=)(2h O =所以，显式Euler 方法是1阶方法，其截断误差的主项是)(22n x y h ''。
7.对隐式Euler 公式推导局部截断误差及其主项,并指出该方法是几阶方法。
解：其局部截断为))(,()()(1111++++--=n n n n n x y x hf x y x y T 对在处作Talor 展开，有)(1+n x y n x )()(2)()()(321h O x y h x y h x y x y n n n n +''+'+=+而且，也在处作Talor 展开，有 ))(,()(111+++='n n n x y x f x y n x )()()()(21h O x y h x y x y n n n +''+'='+所以，因此其局部截断为))(,()()(1111++++--=n n n n n x y x hf x y x y T)()(2)()(32h O x y h x y h x y n n n +''+'+=)()()()(32h O x y h x y h x y n n n +''-'-- )()(232h O x y h n +''-=)(2h O =所以，隐式Euler 方法也是1阶方法，其截断误差的主项是)(22n x y h ''-。
8.对梯形公式推导局部截断误差及其主项,并指出该方法是几阶方法. 解：其局部截断为))](,())(,([2)()(1111+++++--=n n n n n n n x y x f x y x f hx y x y T对在处作Talor 展开，有)(1+n x y n x )()(6)(2)()()(4321h O x y h x y h x y h x y x y n n n n n +'''+''+'+=+而且，))(,()(n n n x y x f x y ='))(,()(111+++='n n n x y x f x y ，对)(1+'n x y 也在处作Talor 展开，有n x )()(2)()()(321h O x y h x y h x y x y n n n n +'''+''+'='+所以，因此其局部截断为))](,())(,([2)()(1111+++++--=n n n n n n n x y x f x y x f hx y x y T)()(6)(2)()(432h O x y h x y h x y h x y n n n n +'''+''+'+=)()(12)(2)(2)(2)(432h O x y h x y h x y h x y h x y n n n n n +'''-''-'-'--)()(1243h O x y h n +'''-=)(3h O =所以，梯形公式是2阶方法，其截断误差的主项是)(123n x y h '''-=。
9.用欧拉法解初值问题，取步长h =0.2.计算过程保留4位小数.⎩⎨⎧1=060≤≤0--='2)().(y x xy y y 解: h =0.2, f (x )=－y －xy 2.首先建立欧拉迭代公式)2,1,0)(4(2.0),(21=-=--=+=+k y x y y hx hy y y x hf y y k k k kk k k k k k k 当k =0，x 1=0.2时，已知x 0=0,y 0=1，有y (0.2)≈y 1=0.2×1(4－0×1)＝0.800 0当k ＝1，x 2=0.4时，已知x 1=0.2, y 1=0.8，有 y (0.4)≈y 2=0.2×0.8×(4－0.2×0.8)＝0.614 4 当k =2，x 3=0.6时，已知x 2=0.4,y 2=0.614 4，有 y (0.6)≈y 3=0.2×0.614 4×(4－0.4×0.4613)＝0.800 010.用欧拉预报－校正公式求解初值问题，取步长h =0.2,计算y (0.2),y (0.4)的近似值，计算过程保留5位小数.l ⎩⎨⎧1=10=++'2)(sin y x y y y 解 步长h =0.2, 此时f (x ,y )=－y －y 2sin x . 欧拉预报－校正公式为：⎪⎩⎪⎨⎧++=+=++++)],(),([2),(1111k k k k k k k k k k y x f y x f hy y y x hf y y 校正值预报值 有迭代公式：⎪⎪⎪⎩⎪⎪⎪⎨⎧+--=--+--+=-=--+=++++++++)sin (1.0)sin 1.09.0()]sin ()sin [(2)sin 2.08.0()sin (121112112121k k k k k k k k k k k k k k k k k k k k k k x y y x y y x y y x y y h y y x y y x y y h y y 校正值预报值当k =0，x 0=1, y 0=1时，x 1=1.2，有631710=11⨯02-80⨯1=20-80=0001.)sin .()sin ..(x y y y715490=21631710+63171010-1⨯1⨯10-90⨯1=≈2121.).sin ..(.)sin ..().(y y 当k =1，x 1=1.2, y 1=0.71549时，x 2=1.4，有476970=21715490⨯02-80⨯715490=20-80=1112.).sin ..(.)sin ..(x y y y).sin ..(.).sin ...(.).(41476970+47697010-21⨯715490⨯10-90⨯715490=≈4122y y =0.5260811.用改进的欧拉法平均公式，取步长h =0.1，求解初值问题⎩⎨⎧=≤≤+='1)0()2.00(y x y x y 计算过程保留4位小数. 解 首先建立迭代格式：⎪⎪⎪⎩⎪⎪⎪⎨⎧+++++=+=++++=+=++=+=++++k k k c p k kk k p k k c k k k k k p y h h hx x h h y y y x h hx h h y y x hf y y h y hx y x hf y y )21(])1([21][21)1(),()1(),(2112121 当k =0时，x 0=0,y 0=1,x 1=0.1,有11.11)21.01.01(]1.01.00)1.01(1.0[2121=⨯+++⨯+⨯+⨯=y当k =1时，x 1=0.1, y 1=1.11, x 2=0.2,有1242.111.121.01.01(]2.01.01.0)1.01(1.0[2122=⨯+++⨯+⨯+⨯=y12.（1）取步长h=0.2,用改进Euler 法求解常微分方程初值问题22,(0)y x y y '=++=0在x=0.6上的解。
（2）对改进Eluer 格式进行误差分析。
解:（1）改进Euler 公式⎪⎩⎪⎨⎧++=+=++++))~,(),((2),(~1111i i i i i i i i i i y x f y x f h y y y x hf y y 分别将x=0.2, 0.4,0.6,代入上式中计算即可！(2)改进欧拉格式11121[]2(,)(,i i i i i i h 2)y y k k k f x y k f x h y hk +⎧=++⎪⎪=⎨⎪=++⎪⎩2341()()()()()()()2!3!i i i i i i h h y x y x h y x hy x y x y x O h +''''''∴=+=++++1211222211(,)()(,)(,)(,)(,1[(,)2(,)(,)]()2!i i i i i i i x i i y i i xxi i xy i i yy i i k f x y y x k f x h y hk f x y hf x y hk f x y h f x y h k f x y h k f x y O h '==''=++=++''''''++++3)又1()(,)(,)()(,)(,)i x i i y i i i x i i y i i y x f x y f x y y x f x y f x y k '''''''=+=+代入11[2i i h 2]y y k k +=++，整理后 31()()(i i )y x y x O h +-=13. （1）取步长为0.1，试用欧拉公式求解常微分方程初值问题⎩⎨⎧=-+='1)0(1y y x y 在x =0.4处的近似值(计算过程保留3位小数)；(2) 试用泰勒展式估计改进欧拉公式的局部截断误差。