【精选】第8章 常微分方程初值问题的数值解法
常微分方程初值问题的数值解法
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
第8章-常微分方程初值问题数值解法
一个数值方法。
h
yi1 yi 2 f (xi , yi ) f (xi1, yi1 ) ( 8.5 )
(8.5)式的右端含有未知的yi+1,它是一个关于 yi+1的函数方程,这类数值方法称为隐式方法。相 反地,欧拉法是关于yi+1的一个直接的计算公式， 这类数值方法称为显式方法。
对于欧拉公式，假定 yi y(xi )，则有
得梯形公式（隐式公式）。
例8.1 用欧拉法解初值问题
y y xy2 (0 x 0.6) y(0) 1
取步长h=0.2 ,计算过程保留4位小数 解: h=0.2, f (x, y) y xy2 欧拉迭代格式
yi1 yi hf (xi , yi ) yi hyi hxi yi2
yn
h( yn
2xn ) yn
1.1yn
0.2xn . yn
计算及结果如下.
clear; y=1, x=0, %初始化
y= 1
x= 0
y = 1.1000 x = 0.1000
for n=1:10
y = 1.1918 x = 0.2000
y=1.1*y-0.2*x/y, x=x+0.1,
方程为
y y0 f (x0 , y0 )( x x0 )
当x x1时,得
y1 y0 f (x0 , y0 )( x1 x0 )
这样就获得了P1点的坐标。
P1
P1 P0
Pi+1 Pn Pi
Pi Pi+1
y=y(x) Pn
x0 x1
xi xi+1 xn
同样, 过点P1(x1,y1),作积分曲线y=y(x)的切线
第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展开式有较多的项达到一致，从而得出较高阶的数 值公式，这就是龙格—库塔法的基本思想。
常微分方程初值问题数值解法
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题，如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题，解析解难以求得或不存在，因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法，以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术，提高计算效率，处理大规模问题。
自适应方法
研究自适应算法，根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法（如欧拉方法）得到一个 初步解，然后使用校正方法（如龙格-库塔方法） 对初步解进行修正，以提高精度。
预估校正方法的优点是精度较高，且计算量相 对较小，适用于各种复杂问题。
步长与误差控制
01
在离散化过程中，步长是一个重要的参数，它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节，它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值，方便快捷，适用范围广。
数值解法的历史与发展
早期发展
早在17世纪，科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展，数值 解法在20世纪取得了重要进展， 如欧拉法、龙格-库塔法等。
当前研究热点
目前，常微分方程初值问题的 数值解法仍有许多研究热点和 挑战，如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。
第8章 常微分方程初值问题数值解法
略去的项称为局部截断误差, 即假定yi y ( xi )时, 近似值yi 1与精确值y ( xi +1 )的差 若局部截断误差R[ y ] O(h p 1 ), 则称该数值方法是p阶方法。
例81 . 用上述三种方法求初值问题 y 0.9 y /(1 2 x), y(0) 1 从x 0到x 1的数值解，取步长h 0.02.
1 yn 1 yn (k1 k2 ) 2 2 xn 改进的Euler法： ) k1 0.1( yn yn 2( xn 0.1) k 0.1( y k ) 2 n 1 yn k1 计算结果如下表所示：
(n 0,1, 2, )
h2 y ''(i ), y ( xi 1 ) y ( xi ) hf ( xi , y ( xi )) 2 h2 y( xi 1) y( xi ) hf ( xi 1, y( xi 1)) y ''(i ) 2 h3 y ( xi 1 ) y ( xi 1 ) 2hf ( xi , y ( xi )) y '''(i ) 3
但不是该公式结构整齐算式求得经典四阶类似以上推导过程可623rk法的稳定性以经典rk法的稳定性相应的误差为于是rk法的稳定性稳定性必须很小才能保证算法限制很大步长的绝对值较大时因此可得稳定性区间为负实数时于是得稳定性区域为由稳定性要求又能保证数值稳定精度提高既可以构造高阶隐式方法因此限制较小即对步长它们均是但从稳定性分析知迭代求解是它们的缺点需要多次法和梯形方法计算隐式格式的向后否则数值不稳定有较大的限制对步长较大时但从稳定性分析知使用方便可以直接计算出法的优点是从22211211约定方法形式为
微分方程初值问题的数值解法
积分法:
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
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解y′=f(x,y),x∈[x0,b]y(x0)=y0.定理1（利普希茨条件）若存在正数L，使得对任意，y1，y2，有|f(x,y1)−f(x,y2)|≤L|(y1−y2)|定理2（解存在性）①若函数f在方区域x∈[a,b],y∈R连续，②函数f关于y 满足利普希茨条件，则对任意x∈[a,b]，常微分方程存在唯一的连续可微数值解.两类问题：①单步法---计算下一个点的值yn+1只需要用到前面一个点的值yn②多步法---计算下一个点的值yn+1需要用到前面l个点的值yl1、欧拉法---下一个点的计算值等于前一个点的计算值加上步长乘以前一个点的函数值•具体过程一些批注：显式欧拉方程指下一步要计算的值，不在迭代方程中；隐式欧拉方程指下一步要计算的值，在迭代方程中。
怎么计算隐式欧拉方程----要借助显示欧拉迭代计算---一般用迭代法-----迭代---将微分方程在区间[xn,xn+1]进行积分，然后函数f进行近似，即可得到迭代方程-----迭代方程收敛性？由函数关于y满足利普希茨条件，可以推出迭代公式收敛。
•局部截断误差：假设前n步误差为0，我们计算第n+1步的误差，将次误差称为局部截断误差，且局部误差为O(hp+1)•p阶精度：由理论证明：若局部误差阶的时间复杂度为O(hp+1)，则整体误差阶为O(hp)我们称公式精度为p。
•显示欧拉法与隐式欧拉法•梯形方法----将显式欧拉迭代方程与隐式欧拉迭代方程做一下加权平均，构造的计算公式.•改进的欧拉方法---思想：因为梯形公式是隐式公式，将显式欧拉公式对下一步的计算值进行预估，用梯形公式对下一步的计算值进行校正.2、龙格-库塔方法思想：根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以前一个点的斜率；而这个斜率用该区间上的多个点的斜率的算数平均来逼近。
注意：怎么计算任意斜率Ki？第i个点的斜率Ki有微分方程可以算出f′=f(xn,yn)所以要算的f(xn,yn)值，由欧拉法即可算出, yn+1=yn+hf′•2阶-龙格-库塔方法----类似改进的欧拉法根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
第8章--常微分方程边值问题的数值解法
第8章 常微分方程边值问题的数值解法8.1 引 言推论 若线性边值问题()()()()()(),,(),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤⎧⎨==⎩ (8.1.2) 满足(1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。
求边值问题的近似解，有三类基本方法：(1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解；(2) 有限元法(finite element method)；(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。
8.2 差分法8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法设二阶线性常微分方程的边值问题为(8.2.1)(8.2.2)()()()(),,(),(),y x q x y x f x a x b y a y b αβ''-=<<⎧⎨==⎩其中(),()q x f x 在[,]a b 上连续，且()0q x ≥.用差分法解微分方程边值问题的过程是：(i) 把求解区间[,]a b 分成若干个等距或不等距的小区间，称之为单元；(ii) 构造逼近微分方程边值问题的差分格式. 构造差分格式的方法有差分法, 积分插值法及变分插值法；本节采用差分法构造差分格式；(iii) 讨论差分解存在的唯一性、收敛性及稳定性；最后求解差分方程. 现在来建立相应于二阶线性常微分方程的边值问题(8.2.1), (8.2.2)的差分方程. ( i ) 把区间[,]I a b =N 等分，即得到区间[,]I a b =的一个网格剖分：011N N a x x x x b -=<<<<=,其中分点(0,1,,)i x a ih i N =+=，并称之为网格节点(grid nodes)；步长b a Nh -=. ( ii ) 将二阶常微分方程(8.2.2)在节点i x 处离散化：在内部节点(1,2,,1)i x i N =-处用数值微分公式2(4)1112()2()()()(),12i i i i i i i i y x y x y x h y x y x x h ξξ+---+''=-<< (8.2.3)代替方程(8.2.2)中()i y x ''，得112()2()()()()()()i i i i i i i y x y x y x q x y x f x R x h +--+-=+,(8.2.4) 其中2(4)()()12i i h R x y ξ=. 当h 充分小时，略去式(8.2.4)中的()i R x ，便得到方程(8.2.1)的近似方程1122i i i i i i y y y q y f h +--+-=,(8.2.5)其中(),()i i i i q q x f f x ==, 11,,i i i y y y +-分别是11(),(),()i i i y x y x y x +-的近似值, 称式(8.2.5)为差分方程(difference equation)，而()i R x 称为差分方程(8.2.5)逼近方程(8.2.2)的截断误差(truncation error). 边界条件(8.7.2)写成0,.N y y αβ==(8.2.6)于是方程(8.2.5), (8.2.6)合在一起就是关于1N +个未知量01,,,N y y y ，以及1N +个方程式的线性方程组：2211212211222111(2),(2),1,2,,1,(2).i i i i i N N N N q h y y h f y q h y y h f i N y q h y h f αβ-+----⎧-++=-⎪-++==-⎨⎪-+=-⎩(8.2.7)这个方程组就称为逼近边值问题(8.2.1), (8.2.2)的差分方程组(system of difference equations)或差分格式(difference scheme)，写成矩阵形式2211122222223332222222111(2)11(2)11(2)11(2)11(2)N N N N N N y q h h f y q h h f y q h h f y q h h f y q h h f αβ------⎡⎤⎡⎤-+-⎡⎤⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦. (8.2.8)用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.7)或(8.2.8), 其解01,,,N y y y 称为边值问题(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)，我们自然关心它是否有唯一解；此外，当网格无限加密，或当0h →时，差分解i y 是否收敛到微分方程的解()i y x . 为此介绍下列极值原理：定理8.2.1 (极值原理) 设01,,,N y y y 是给定的一组不全相等的数，设1122(),0,1,2,,i i i i i i i y y y l y q y q i N h +--+=-≥=.(8.2.9)(1) 若()0,1,2,,i l y i N ≥=, 则{}0Ni i y =中非负的最大值只能是0y 或N y ； (2) 若()0,1,2,,i l y i N ≤=, 则{}0Ni i y =中非正的最小值只能是0y 或N y .证 只证(1) ()0i l y ≥的情形，而(2) ()0i l y ≤的情形可类似证明. 用反证法. 记{}0max i i NM y ≤≤=，假设0M ≥, 且在121,,,N y y y -中达到. 因为i y 不全相等，所以总可以找到某个00(11)i i N ≤≤-，使0i y M =，而01i y -和01i y +中至少有一个是小于M 的. 此时0000000011222()2.i i i i i i i i y y y l y q y h M M M q M q M h +--+=--+<-=-因为00,0i q M ≥≥，所以0()0i l y <, 这与假设矛盾，故M 只能是0y 或N y . 证毕！推论 差分方程组(8.2.7)或(8.2.8)的解存在且唯一. 证明 只要证明齐次方程组11202()0,0,1,2,,1,0,0i i i i i i i N y y y l y q y q i N h y y +--+⎧=-=≥=-⎪⎨⎪==⎩ (8.2.10)只有零解就可以了. 由定理8.7.1知，上述齐次方程组的解01,,,N y y y 的非负的最大值和非正的最小值只能是0y 或N y . 而00,0N y y ==，于是0,1,2,,.i y i N == 证毕！利用定理8.2.1还可以证明差分解的收敛性及误差估计. 这里只给出结果： 定理8.2.2 设i y 是差分方程组(8.2.7)的解，而()i y x 是边值问题(8.2.1), (8.2.2)的解()y x 在i x 上的值，其中0,1,,i N =. 则有224()(),96i i i M h y x y b a ε=-≤-(8.2.11)其中(4)4max ()a x bM yx ≤≤=.显然当0h →时，()0i i i y x y ε=-→. 这表明当0h →时，差分方程组(8.2.7)或(8.2.8)的解收敛到原边值问题(8.7.1), (8.7.2)的解.例8.2.1 取步长0.1h =，用差分法解边值问题43,01,(0)(1)0,y y x x y y ''-=≤≤⎧⎨==⎩并将结果与精确解()()2222()3434x xy x e e ee x --=---进行比较.解 因为110N h ==，()4,()3q x f x x ==, 由式(8.2.7)得差分格式221222112289(240.1)30.10.1,(240.1)30.1,2,3,,8,(240.1)30.10.9,i i i i y y y y y x i y y -+⎧-+⨯+=⨯⨯⎪-+⨯+=⨯=⎨⎪-+⨯=⨯⨯⎩0100y y ==, 00.1,1,2,,9i x ih i i =+==, 其结果列于表8.2.1.从表8.2.1可以看出, 差分方法的计算结果的精度还是比较高的. 若要得到更精确的数值解，可用缩小步长h 的方法来实现.8.2.2 一般二阶线性常微分方程边值问题的差分法对一般的二阶微分方程边值问题1212()()()()()(),,()(),()(),y x p x y x q x y x f x a x b y a y a y b y b αααβββ'''++=<<⎧⎪'+=⎨⎪'+=⎩ (8.2.12) 假定其解存在唯一.为求解的近似值，类似于前面的做法，( i ) 把区间[,]I a b =N 等分，即得到区间[,]I a b =的一个网格剖分：011N N a x x x x b -=<<<<=,其中分点(0,1,,)i x a ih i N =+=，步长b a Nh -=. ( ii ) 对式(8.2.12)中的二阶导数仍用数值微分公式2(4)1112()2()()()(),12i i i i i i i iy x y x y x h y x y x x h ξξ+---+''=-<<代替，而对一阶导数，为了保证略去的逼近误差为2()O h ，则用3点数值微分公式；另外为了保证内插，在2个端点所用的3点数值微分公式与内网格点所用的公式不同，即21112012000022212()()()(),,1,2,,1,263()4()()()(),,23()4()3()()(),.23i i i i i i i N N N N N N N N y x y x h y x y x x i N h y x y x y x h y x y x x h y x y x y x h y x y x x h ξξξξξξ+-----⎧-''''=-<<=-⎪⎪-+-⎪''''=+<<⎨⎪⎪-+''''=+<<⎪⎩(8.2.13) 略去误差，并用()i y x 的近似值i y 代替()i y x ，(),(),()i i i i i i p p x q q x f f x ===，便得到差分方程组1111221001221211(2)(),1,2,,1,2(34),2(43),2i i i i i i i i i N N N N p y y y y y q y f i N h hy y y y h y y y y h αααβββ-++---⎧-++-+==-⎪⎪⎪+-+-=⎨⎪⎪+-+=⎪⎩(8.2.14)其中(),(),(),1,2,,1i i i i i i q q x p p x f f x i N ====-, i y 是()i y x 的近似值. 整理得12021222211222121(23)42,(2)2(2)(2)2,1,2,,1,4(32)2.i i i i i i i N N N h y y y h hp y h q y hp y h f i N y y h y h αααααβββββ-+---+-=⎧⎪---++==-⎨⎪-++=⎩ (8.2.15)解差分方程组(8.2.15)，便得边值问题(8.2.12)的差分解01,,,N y y y .特别地, 若12121,0,1,0ααββ====，则式(8.2.12)中的边界条件是第一类边值条件：(),();y a y b αβ==此时方程组(7.7.16)为221112112211221211112(2)(2)2(2),(2)2(2)(2)2,2,3,,2,(2)2(2)2(2).i i i i i i i N N N N N N h q y hp y h f hp hp y h q y hp y h f i N hp y h q y h f hp αβ-+------⎧--++=--⎪---++==-⎨⎪---=-+⎩(8.2.16) 方程组(8.2.16)是三对角方程组，用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.16)，便得边值问题(8.2.12)的差分解01,,,N y y y .( iii ) 讨论差分方程组(8.2.16)的解是否收敛到微分方程的解，估计误差. 这里就不再详细介绍.例8.2.2取步长/16h π=，用差分法求下列边值问题的近似解，并将结果与精确解进行比较.精确解是1()(sin 3cos )10y x x x =-+. 解 因为(20)8N h π=-=，()1,()2,()cos p x q x f x x =-=-=, 由式(8.2.17)得差分格式()()()()()()()()()()()()()2122211222122216(2)216(1)216cos 16216(1)(0.3),216(1)2216(2)216(1)216cos 16,2,3,,6,216(1)2216(2)216cos 7i i i N N y yy y y i i y y πππππππππππππ-+--⎡⎤--⨯-++⨯-⎡⎤⎣⎦⎣⎦=--⨯-⨯-⎡⎤⎣⎦⎡⎤-⨯---⨯-++⨯-⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦==⎡⎤-⨯---⨯-⎡⎤⎣⎦⎣⎦=()()16216(1)(0.1),ππ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪-+⨯-⨯-⎡⎤⎣⎦⎩080.3,0.1y y =-=-, 00.1,1,2,,9i x ih i i =+==, 其结果列于表8.2.2.8.3 有限元法有限元法(finite element method)是求解微分方程定解问题的有效方法之一，它特别适用在几何、物理上比较复杂的问题. 有限元法首先成功地应用于结构力学和固体力学，以后又应用于流体力学、物理学和其他工程科学. 为简明起见，本节以线性两点边值问题为例介绍有限元法.考虑线性两点边值问题()(8.3.1)(8.3.2)()()()()(),,(),(),Ly p x y x q x y x f x a x b y a y b αβ⎧''⎪=-+=≤≤⎨==⎪⎩其中1()0,()0,C [,]p x q x p a b >≥∈, ,C[,]q f a b ∈.此微分方程描述了长度为b a -的可变交叉截面(表示为()q x )的横梁在应力()p x 和()f x 下的偏差()y x .8.3.1 等价性定理记{}221C [,]()C [,],(),()a b y y y x a b y a y b αβ==∈==, 引进积分()22()()[()]()()2()()d baI y p x y x q x y x f x y x x '=+-⎰. (8.3.3)任取21()C [,]y y x a b =∈，就有一个积分值()I y 与之对应，因此()I y 是一个泛函(functional)，即函数的函数. 因为这里是,y y '的二次函数，因此称()I y 为二次泛函.对泛函(8.3.3)有如下变分问题(variation problem)：求函数21C [,]y a b ∈，使得对任意21C [,]y a b ∈, 均有()()I y I y ≥, (8.3.4)即()I y 在y 处达到极小, 并称y 为变分问题(8.3.4)的解.可以证明：定理8.3.1(等价性定理) y 是边值问题(8.3.1), (8.3.2)的解的充分必要条件是y 使泛函()I y 在21C [,]a b 上达到极小，即y 是变分问题(8.3.4)在21C [,]a b 上的解. 证 (充分性) 设21C [,]y a b ∈是变分问题()I y 的解；即y 使泛函()I y 在21C [,]a b 上达到极小，证明y 必是边值问题(8.3.1), (8.3.2)的解.设()x η是2C [,]a b 任意一个满足()()0a b ηη==的函数，则函数21()()()C [,]y x y x x a b αη=+∈，其中α为参数. 因为y 使得()I y 达到极小，所以()()I y I y αη+≥，即积分()22()()[()()]()[()()]2()[()()]baI y p x y x x q x y x x f x y x x dxαηαηαηαη''+=+++-+⎰作为α的函数，在0α=处取极小值()I y ，故d()0d I y ααηα=+=. (8.3.5) 计算上式，得()()()()()022(8.d()d d ()[()()]()[()()]2()[()()]d d 2()[()()]()2()[()()]()2()()d 2()()()()()()()()d .bab abaI y p x y x x q x y x x f x y x x x p x y x x x q x y x x x f x x x p x y x x q x y x x f x x x ααααηααηαηαηααηηαηηηηηη===+''=+++-+'''=+++-''=+-⎰⎰⎰3.6)利用分部积分法计算积分[][]()()()d ()()d ()()()()()()()d ()()()d ,bbaab ba abap x y x x x p x y x x p x y x x x p x y x x x p x y x x ηηηηη'''='''=-''=-⎰⎰⎰⎰代入式(8.3.6)，得()0(8.3.7)d()2()()()()()()d 0.d b a I y p x y x q x y x f x x x ααηηα'=⎡⎤⎣⎦'+=-+-=⎰因为()x η是任意函数，所以必有()()()()()()0p x y x q x y x f x ''-+-≡. (8.3.8)否则，若在[,]a b 上某点0x 处有()00000()()()()()0p x y x q x y x f x ''-+-≠，不妨设()00000()()()()()0p x y x q x y x f x ''-+->，则由函数的连续性知，在包含0x 的某一区间00[,]a b 上有()()()()()()0p x y x q x y x f x ''-+->.作002200000,[,]\[,],()()(),.x a b a b x x a x b a x b η∈⎧⎪=⎨--≤≤⎪⎩ 显然2()C [,]x a b η∈，且()()0a b ηη==，但()()()()()()()d ba p x y x q x y x f x x x η⎡⎤''-+-⎢⎥⎣⎦⎰ ()00()()()()()()d 0b a p x y x q x y x f x x x η⎡⎤''=-+->⎢⎥⎣⎦⎰，这与式(8.3.7)矛盾. 于是式(8.3.8)成立，即变分问题(8.3.4)的解y 满足微分方程(8.3.1), 且(),()y a y b αβ==故它是边值问题(8.3.1), (8.3.2)的解.(必要性) 设()y y x =是边值问题(8.3.1), (8.3.2)的解，证明y 是变分问题(8.3.4)的解；即：y 使泛函()I y 在21C [,]a b 上达到极小.因为()y y x =满足方程(8.3.1)，所以()()()()()()0p x y x q x y x f x ''-+≡.设任意21()C [,]y y x a b =∈，则函数()()()x y x y x η=-满足条件()()0a b ηη==，且2()C [,]x a b η∈. 于是()()[]()222222()()()()()[()()]()[()()]2()[()()]d ()[()]()[()]2()()d 2()()()()()()()()d ()[()]()[()]d baba baaI y I y I y I y p x y x x q x y x x f x y x x x p x y x q x y x f x y x xp x y x x q x y x x f x x x p x x q x x xηηηηηηηηη-=+-''=+++-+'-+-''=+-++⎰⎰⎰()()()22222()()()()()()d ()[()]()[()]d ()[()]()[()]d .bb ba a bap x y x q x y x f x x x p x x q x x xp x x q x x x ηηηηη⎡⎤'''=--+++⎢⎥⎣⎦'=+⎰⎰⎰⎰因为()0,()0p x q x >≥，所以当()0x η≠时，()22()[()]()[()]d 0bap x x q x x x ηη'+>⎰, 即()()0I y I y ->.只有当()0x η≡时，()()0I y I y -=. 这说明y 使泛函()I y 在21C [,]a b 上达到极小. 证毕！定理8.3.2 边值问题(8.3.1), (8.3.2)存在唯一解.证明 用反证法. 若12(),()y x y x 都是边值问题(8.3.1), (8.3.2)的解，且不相等，则由定理8.3.1知，它们都使泛函()I y 在21C [,]a b 上达到极小，因而12()()I y I y > 且 21()()I y I y >，矛盾！因此边值问题(8.3.1), (8.3.2)的解是唯一的.由边值问题解的唯一性，不难推出边值问题(8.3.1), (8.3.2)解的存在性(留给读者自行推导).8.3.2 有限元法等价性定理说明，边值问题(8.3.1), (8.3.2)的解可化为变分问题(8.3.4)的求解问题. 有限元法就是求变分问题近似解的一种有效方法. 下面给出其解题过程：第1步 对求解区间进行网格剖分01,i n a x x x x b =<<<<<=区间1[,]i i i I x x -=称为单元，长度1(1,2,,)i i i h x x i n -=-=称为步长，1max i i nh h ≤≤=. 若(1,2,,)i h h i n ==，则称上述网格剖分为均匀剖分.给定剖分后，泛函(8.3.3)可以写成()22()()[()]()()2()()d baI y p x y x q x y x f x y x x '=+-⎰()12211()[()]()()2()()d i i nnx i x i i p x y x q x y x f x y x xS -=='=+-∑∑⎰记. (8.3.9)第2步 构造试探函数空间。
计算方法课件第八章常微分方程初值问题的数值解法
整体截断误差与局部截断误差的关系
定理：如果f(x，y)满足李普希兹（Lipschitz）条件
f(x ,y 1 )f(x ,y 2) L y 1y 2
且局部截断误差有界：
|R n|1 2h2M 2
(n1,2, )
则Euler法的整体截断误差n满足估计式：
ne(ba)L 0h 2L M 2(e(ba)L1)
分光滑。初值问题的解析解(理论解）用 y(x表n ) 示, 数值解法的精确解用 y表n 示。
常微分方程数值解法一般分为：
（1）一步法：在计算y n 1 时，只用到x n 1 ,x n和 y，n 即前一步的值。
（2）多步法：计算 y n 1 时，除用到 x n 1 ,x n 和 y n 以外，还要用 x n p 和 y n p (p1 ,2 k;k0) ,即前
其中L为李普希兹常数，b-a为求解区间长度，
M2 mayx(x) 。 axb
证明参见教材。
Remark：该定理表明，整体截断误差比局部截 断误差低一阶。对其它方法，也有类似的结论。
收敛性与稳定性
收敛性定义：如果某一数值方法对于任意固定的
xn=x0+nh，当h0(同时n )时有yn y(xn)，
则称该方法收敛。 稳定性定义 定义 用一个数值方法，求解微分方程初值问 题时，对给定的步长h>0，若在计算 y n 时引入 误差 （n 也称扰动），但由此引起计算后面的 ynk(k1,2, )时的误差按绝对值均不增加，则 称这个数值方法是稳定的。
一般的显式rk方法可以写成型钢截面只需少量加工即可用作构件省工省时成本低但型钢截面受型钢种类及型钢号限制难于完全与受力所需的面积相对应用料较多其中为常数选取这些常数的原则是要求第一式的右端在处泰勒展开后按h型钢截面只需少量加工即可用作构件省工省时成本低但型钢截面受型钢种类及型钢号限制难于完全与受力所需的面积相对应用料较多上述公式叫做n级的rungekutta方法其局部截断误差为显然euler法是一级一阶rk方法
常微分方程数值解法
第八章 常微分方程的数值解法一．内容要点考虑一阶常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。
在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。
用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。
(一)常微分方程处置问题解得存在唯一性定理对于常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy如果：(1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。
(2) ),(y x f 对于y 满足利普希茨条件，即2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。
定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。
收敛性定理：若一步方法满足： (1)是p 解的.(2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件.(3) 初始值y 0是精确的。
则),()()(p h O x y kh y =-kh =x -x 0,也就是有0x y y lim k x x kh 0h 0=--=→)((一)、主要算法 1．局部截断误差局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~+k y 的误差y (x k+1)- 1~+k y 称为局部截断误差。
常微分方程初值问题的数值解法
1
1、基本概念： 基本概念：
•一阶常微分方程初值问题是： 一阶常微分方程初值问题是： 一阶常微分方程初值问题是 y ′ =f(x,y) (1.1) y(x0)=y0 (1.2) 其中f是已知的xy平面上某个区域D上连续函数，式（1.1）是微分 方程，有无穷多解，式（1.2）是确定解的初始条件。 如一元函数y(x)对一切a≤x≤b 满足 （1） （x,y(x))∈D （2） y(x0)=y0 （3） y′存在，且y′(x)=f(x,y(x)) 则称y(x)是初值问题（1.1）、（1.2）在[a,b]上的解。 定义1.1 如果存在正常数L>0,使得对一切x∈[a,b]及y1 、y2 ∈[c,d] 有 |f(x,y1)-f(x,y2)|≤ L|y1- y2| (1.3) 则称f(x,y)满足对y的lipschitz(李普希滋）条件，其中L称为 lipschitz常数。
~ = y + hf ( x , y ) y k +1 k k k y k +1
2011-12-12
h y = y k + ( f ( x k , y k ) + f ( x k +1 , ~k +1 )) 2
12
Euler公式的几何意义 公式的几何意义
Y=y(x)
a
b
2011-12-12
13
例题：用Euler的数值解(取步长h=0.1计算到y3)：
2011-12-12 4
求初值问题，是给出它的解在某些节点数值的近似值，这称为 数值离散方法，若要求出在[a,b} 上的离散解，引入点列{xK} ， 其中 xk= xk-1 + hk-1 ，k=1，2，3….. xk称为节点， hk 称为步长。 通常，步长h不变，取为等距步长 hk =h= （b-a）/N ，N为等份 区间[a,b]分割数,节点变为等距节点，此时有 kh, xk= xk-1 + h= x0 + kh 记式（1.1），（1.2）的准确解y(x)在节点xk 之值为y(xk) ，而 记y(xk)的近似值为yk ，又记fk=f（xk ， yk ），
第八章 常微分方程初值问题的解法
第八章常微分方程初值问题的解法在科学与工程问题中，常微分方程描述物理量的变化规律，应用非常广泛. 本章介绍最基本的常微分方程初值问题的解法，主要针对单个常微分方程，也讨论常微分方程组的有关技术.8.1引言本节介绍常微分方程、以及初值问题的基本概念，并对常微分方程初值问题的敏感性进行分析.8.1.1 问题分类与可解性很多科学与工程问题在数学上都用微分方程来描述，比如，天体运动的轨迹、机器人控制、化学反应过程的描述和控制、以及电路瞬态过程分析，等等. 这些问题中要求解随时间变化的物理量，即未知函数y(t)，t表示时间，而微分方程描述了未知函数与它的一阶或高阶导数之间的关系. 由于未知函数是单变量函数，这种微分方程被称为常微分方程（ordinary differential equation, ODE），它具有如下的一般形式①：g(t,y,y′,⋯,y(k))=0 ,(8.1) 其中函数g: ℝk+2→ℝ. 类似地，如果待求的物理量为多元函数，则由它及其偏导函数构成的微分方程称为偏微分方程（partial differential equation, PDE）. 偏微分方程的数值解法超出了本书的范围，但其基础是常微分方程的解法.在实际问题中，往往有多个物理量相互关联，它们构成的一组常微分方程决定了整个系统的变化规律. 我们先针对单个常微分方程的问题介绍一些基本概念和求解方法，然后在第8.5节讨论常微分方程组的有关问题.如公式(8.1)，若常微分方程包含未知函数的最高阶导数为y(k)，则称之为k阶常微分方程. 大多数情况下，可将常微分方程(8.1)写成如下的等价形式：y(k)=f(t,y,y′,⋯,y(k−1)) ,(8.2) 其中函数f: ℝk+1→ℝ. 这种等号左边为未知函数的最高阶导数y(k)的方程称为显式常微分方程，对应的形如(8.1)式的方程称为隐式常微分方程.通过简单的变量代换可将一般的k阶常微分方程转化为一阶常微分方程组. 例如对于方程(8.2)，设u1(t)=y(t),u2(t)=y′(t),⋯,u k(t)=y(k−1), 则得到等价的一阶显式常微分方程组为：{u1′=u2u2′=u3⋯u k′=f(t,u1,u2,⋯,u k).(8.3)本书仅讨论显式常微分方程，并且不失一般性，只需考虑一阶常微分方程或方程组.例8.1 (一阶显式常微分方程)：试用微积分知识求解如下一阶常微分方程：y′=y .[解] 采用分离变量法进行推导：①为了表达式简洁，在常微分方程中一般省略函数的自变量，即将y(t)简记为y，y′(t)简记为y′，等等.dy dt =y ⟹ dy y=dt , 对两边积分，得到原方程的解为：y (t )=c ∙e t ,其中c 为任意常数.从例8.1看出，仅根据常微分方程一般无法得到唯一的解. 要确定唯一解，还需在一些自变量点上给出未知函数的值，称为边界条件. 一种边界条件设置方法是给出t =t 0时未知函数的值:y (t 0)=y 0 .在合理的假定下，从t 0时刻对应的初始状态y 0开始，常微分方程决定了未知函数在t >t 0时的变化情况，也就是说这个边界条件可以确定常微分方程的唯一解（见定理8.1）. 相应地，称y (t 0)=y 0为初始条件，而带初始条件的常微分方程问题：{y ′=f (t,y ),t ≥t 0y (t 0)=y 0 . (8.4)为初值问题（initial value problem, IVP ）.定理8.1：若函数f (t,y )关于y 满足李普希兹（Lipschitz ）条件，即存在常数L >0，使得对任意t ≥t 0，任意的y 与y ̂，有：|f (t,y )−f(t,y ̂)|≤L |y −y ̂| ,(8.5) 则常微分方程初值问题(8.4)存在唯一的解.一般情况下，定理8.1的条件总是满足的，因此常微分方程初值问题的解总是唯一存在的. 为了更清楚地理解这一点，考虑f (t,y )的偏导数ðf ðy 存在，则它在求解区域内可推出李普希兹条件(8.5)，因为f (t,y )−f (t,y ̂)=ðf ðy (t,ξ)∙(y −y ̂) , 其中ξ为介于y 和y ̂之间的某个值. 设L 为|ðf ðy (t,ξ)|的上界，(8.5)式即得以满足.对公式(8.4)中的一阶常微分方程还可进一步分类. 若f (t,y )是关于y 的线性函数，f (t,y )=a (t )y +b (t ) ,(8.6) 其中a (t ),b (t )表示自变量为t 的两个一元函数，则对应的常微分方程为线性常微分方程，若b (t )≡0, 则为线性齐次常微分方程. 例8.1中的方程属于线性、齐次、常系数微分方程，这里的“常系数”是强调a (t )为常数函数.8.1.2 问题的敏感性对常微分方程初值问题，可分析它的敏感性，即考虑初值发生扰动对结果的影响. 注意这里的结果（解）是一个函数，而不是一个或多个值. 由于实际应用的需要，分析常微分方程初值问题的敏感性时主要关心t →∞时y (t )受影响的情况，并给出有关的定义. 此外，考虑到常微分方程的求解总与数值算法交织在一起、以及历史的原因，一般用“稳定”、“不稳定”等词汇说明问题的敏感性.定义8.1：对于常微分方程初值问题(8.4)，考虑初值y 0的扰动使问题的解y (t )发生偏差的情形. 若t →∞时y (t )的偏差被控制在有界范围内，则称该初值问题是稳定的（stable ），否则该初值问题是不稳定的（unstable ）. 特别地，若t →∞时y (t )的偏差收敛到零，则称该初值问题是渐进稳定的（asymptotically stable ）.关于定义8.1，说明两点：● 渐进稳定是比稳定更强的结论，若一个问题是渐进稳定的，它必然是稳定的. ● 对于不稳定的常微分方程初值问题，初始数据的扰动将使t →∞时的结果误差无穷大. 因此为了保证数值求解的有效性，常微分方程初值问题具有稳定性是非常重要的.例8.2 (初值问题的稳定性): 考察如下“模型问题”的稳定性：{y ′=λy,t ≥t 0y (t 0)=y 0 . (8.7)[解] 易知此常微分方程的准确解为：y (t )=y 0e λ(t−t 0). 假设初值经过扰动后变为y 0+Δy 0，对应的扰动后解为y ̂(t )=(y 0+Δy 0)e λ(t−t 0),所以扰动带来的误差为Δy (t )=Δy 0e λ(t−t 0) .根据定义8.1，需考虑t →∞时Δy (t )的值，它取决于λ. 易知，若λ≤0，则原问题是稳定的，若λ>0，原问题不稳定. 而且当λ<0时，原问题渐进稳定.图8-1分三种情况显示了初值扰动对问题(8.7)的解的影响，从中可以看出不稳定、稳定、渐进稳定的不同含义.对例8.2中的模型问题，若考虑参数λ为一般的复数，则问题的稳定性取决于λ的实部，若Re(λ)≤0, 则问题是稳定的，否则不稳定. 例8.2的结论还可推广到线性、常系数常微分方程，即根据f (t,y )中y 的系数可确定初值问题的稳定性. 对于一般的线性常微分方程(8.6)，由于方程中y 的系数为关于t 的函数，仅能分析t 取某个值时的局部稳定性.例8.3 (局部稳定性): 考察如下常微分方程初值问题的稳定性：{y ′=−10ty,t ≥0y (0)=1 . (8.8)[解] 此常微分方程为线性常微分方程，其中y 的系数为a (t )=−10t . 当t ≥0时，a (t )≤0，在定义域内每个时间点上该问题都是局部稳定的.事实上，方程(8.8)的解析为y (t )=e −5t 2，初值扰动Δy 0造成的结果误差为Δy (t )=Δy 0e −5t 2. 这说明初值问题(8.8)是稳定的.对于更一般的一阶常微分方程(8.4)，由于其中f (t,y )可能是非线性函数，分析它的稳定性非常复杂. 一种方法是通过泰勒展开用一个线性常微分方程来近似它，再利用线性常微分方程稳定性分析的结论了解它的局部稳定性. 具体的说，在某个解函数y ∗(t)附近用一阶泰勒展开近似f (t,y )，f (t,y )≈f (t,y ∗)+ðf ðy(t,y ∗)∙(y −y ∗) 则原微分方程被局部近似为（用符号z 代替y ）: 图8-1 (a) λ>0对应的不稳定问题, (b) λ=0对应的稳定问题, (c) λ<0对应的渐进稳定问题. (a) (b) (c)z′=ðfðy(t,y∗)∙(z−y∗)+f(t,y∗)这是关于未知函数z(t)的一阶线性常微分方程，可分析t取某个值时的局部稳定性. 因此，对于具体的y∗(t)和t的取值，常微分方程初值问题(8.4)的局部稳定性取决于ðfðy(t,y∗)的实部的正负号. 应注意的是，这样得到的关于稳定性的结论只是局部有效的.实际遇到的大多数常微分方程初值问题都是稳定的，因此在后面讨论数值解法时这常常是默认的条件.8.2简单的数值解法与有关概念大多数常微分方程都无法解析求解（尤其是常微分方程组），只能得到解的数值近似. 数值解与解析解有很大差别，它是解函数在离散点集上近似值的列表，因此求解常微分方程的数值方法也叫离散变量法. 本节先介绍最简单的常微分方程初值问题解法——欧拉法（Euler method），然后给出数值解法的稳定性和准确度的概念，最后介绍两种隐格式解法.8.2.1 欧拉法数值求解常微分方程初值问题，一般都是“步进式”的计算过程，即从t0开始依次算出离散自变量点上的函数近似值. 这些离散自变量点和对应的函数近似值记为：t0<t1<⋯<t n<t n+1<⋯y 0,y1,⋯y n,y n+1,⋯其中y0是根据初值条件已知的. 相邻自变量点的间距为 n=t n+1−t n, 称为步长.数值解法通常使用形如y n+1=G(y n+1,y n,y n−1,…,y n−k)(8.9) 的计算公式，其中G表示某个多元函数. 公式(8.9)是若干个相邻时间点上函数近似值满足的关系式，利用它以及较早时间点上函数近似值可算出y n+1. 若公式(8.9)中k=0，则对应的解法称为单步法（single-step method），其计算公式为：y n+1=G(y n+1,y n) .(8.10) 否则，称为多步法（multiple-step method）. 另一方面，若函数G与y n+1无关，即:y n+1=G(y n,y n−1,…,y n−k),则称为显格式方法（explicit method），否则称为隐格式方法（implicit method）. 显然，显格式方法的计算较简单，只需将已得到的函数近似值代入等号右边，则可算出y n+1.欧拉法是一种显格式单步法，对初值问题(8.4)其计算公式为：y n+1=y n+ n f(t n,y n) , n=0,1,2,⋯.(8.11) 它可根据数值微分的向前差分公式（第7.7节）导出. 由于y′=f(t,y)，则y′(t n)=f(t n,y(t n))≈y(t n+1)−y(t n)n,得到近似公式y(t n+1)≈y(t n)+ n f(t n,y(t n))，将其中的函数值换为数值近似值，则得到欧拉法的递推计算公式(8.11). 还可以从数值积分的角度进行推导，由于y(t n+1)=y(t n)+∫y′(s)dst n+1t n =y(t n)+∫f(s,y(s))dst n+1t n,用左矩形公式近似计算其中的积分（矩形的高为s=t n时被积函数值），则有y(t n+1)≈y(t n)+ n f(t n,y(t n)) ,将其中的函数值换为数值近似值，便得到欧拉法的计算公式.例8.4 (欧拉法)：用欧拉法求解初值问题{y ′=t −y +1y (0)=1. 求t =0.5时y (t )的值，计算中将步长分别固定为0.1和0.05.[解] 在本题中，f (t,y )=t −y +1, t 0=0, y 0=1, 则欧拉法计算公式为：y n+1=y n + (t n −y n +1) , n =0,1,2,⋯当步长h=0.1时，计算公式为y n+1=0.9y n +0.1t n +0.1; 当步长h=0.05时，计算公式为y n+1=0.95y n +0.05t n +0.05. 两种情况的计算结果列于表8-1中，同时也给出了准确解y (t )=t +e −t 的结果.表8-1 欧拉法计算例8.4的结果 h=0.1h=0.05 t ny n y (t n ) t n y n t n y n 0.11.000000 1.004837 0.05 1.000000 0.3 1.035092 0.21.010000 1.018731 0.1 1.002500 0.35 1.048337 0.31.029000 1.040818 0.15 1.007375 0.4 1.063420 0.41.056100 1.070320 0.2 1.014506 0.45 1.080249 0.5 1.090490 1.106531 0.25 1.023781 0.5 1.098737 从计算结果可以看出，步长取0.05时，计算的误差较小.在常微分方程初值问题的数值求解过程中，步长 n ,(n =0,1,2,⋯)的设置对计算的准确性和计算量都有影响. 一般地，步长越小计算结果越准确，但计算步数也越多（对于固定的计算区间右端点），因此总计算量就越大. 在实际的数值求解过程中，如何设置合适的步长达到准确度与效率的最佳平衡是很重要的一个问题.8.2.2数值解法的稳定性与准确度在使用数值方法求解初值问题时，还应考虑数值方法的稳定性. 实际的计算过程中都存在误差，若某一步的解函数近似值y n 存在误差，在后续递推计算过程中，它会如何传播呢？会不会恶性增长，以至于“淹没”准确解？通过数值方法的稳定性分析可以回答这些问题. 首先给出稳定性的定义.定义8.2：采用某个数值方法求解常微分方程初值问题(8.4)，若在节点t n 上的函数近似值存在扰动δn ,由它引起的后续各节点上的误差δm (m >n )均不超过δn ，即|δm |≤|δn |,(m >n)，则称该方法是稳定的.在大多数实际问题中，截断误差是常微分方程数值求解中的主要计算误差，因此我们忽略舍入误差. 此外，仅考虑稳定的常微分方程初值问题.考虑单步法的稳定性，需要分析扰动δn 对y n+1的影响，推导δn+1与δn 的关系式. 以欧拉法为例，先考虑模型问题(8.7)，并且设Re(λ)≤0. 此时欧拉法的计算公式为②：y n+1=y n + λy n =(1+ λ)y n ,由y n 上的扰动δn 引起y n+1的误差为：δn+1=(1+ λ)δn ,要使δn+1的大小不超过δn ，则要求|1+ λ|≤1 . (8.12)② 对于稳定性分析以及后面的一些场合，由于只考虑一步的计算，将步长 n 记为 .。
常微分方程初值问题的数值解法
--------学习小结一、本章学习体会通过本章的学习，我了解了常微分方程初值问题的计算方法，对于解决那些很难求解出解析表达式的，甚至有解析表达式但是解不出具体的值的常微分方程非常有用。
在这一章里求解常微分方程的基本思想是将初值问题进行离散化，然后进行迭代求解。
在这里将初值问题离散化的方法有三种，分别是差商代替导数的方法、Taylor 级数法和数值积分法。
常微分方程初值问题的数值解法的分类有显示方法和隐式方法，或者可以分为单步法和多步法。
在这里单步法是指计算第n+1个y 的值时,只用到前一步的值，而多步法则是指计算第n+1个y 的值时，用到了前几步的值。
通过对本章的学习，已经能熟练掌握如何用Taylor 级数法去求解单步法中各方法的公式和截断误差，但是对线性多步法的求解理解不怎么透切，特别是计算过程较复杂的推理。
在本章的学习过程中还遇到不少问题，比如本章知识点多，公式多，在做题时容易混淆，其次对几种R-K 公式的理解不够透彻，处理一个实际问题时，不知道选取哪一种公式，通过课本里面几种方法的计算比较得知其误差并不一样，，这个还需要自己在往后的实际应用中多多实践留意并总结。
二、本章知识梳理常微分方程初值问题的数值解法一般概念步长h ，取节点0,(0,1,...,)n t t nh n M =+=，且M t T ≤，则初值问题000'(,),()y f t y t t T y t y =≤≤⎧⎨=⎩的数值解法的一般形式是1(,,,...,,)0,(0,1,...,)n n n n k F t y y y h n M k ++==-显示单步法7.2.1 显示单步法的一般形式1(,,),(0,1,...,1)n n n n y y h t y h n M ϕ+=+=-定理7.2.1 设增量函数(,,)n n t y h ϕ在区域00{(,,)|,||,0}D t y h t t T y h h =≤≤<∞≤≤内对变量y 满足Lipschitz 条件，即存在常数K ，使对D 内任何两点1(,,)t u h 和2(,,)t u h ，不等式1212|(,,)(,,)|||t u h t u h K u u ϕϕ-≤-成立，那么，若单步法的局部截断误差1n R +与1(1)p h p +≥同阶，即11()p n R O h ++=，则单步法的整体截断误差1n ε+与p h 同阶，即1()p n O h ε+=。
计算方法第八章
y y
( k +1) n +1
(0) n +1
= yn + hf ( xn , yn )
(k ) n +1
= yn + f ( xn , yn ) ⋅ h / 2 + f ( xn +1 , y ) ⋅ h / 2
收敛的条件： 收敛的条件：
f (x, y) = −y yn+1 = yn−1 − 2hyn
y1 = (1− h) y0 , y0 =1
计 算 结 果
0.02 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0.9802 0.90484 0.81873 0.74082 0.67032 0.60653 0.54881 0.49659 0.44933 0.40657
欧拉中点公式
yn+1 = yn−1 + 2hf (xn , yn ) , y(x0 ) = y0
利用中点公式求解微分方程时，有一个问题，就是计算时需 利用中点公式求解微分方程时，有一个问题， 要两个迭代初值！ 要两个迭代初值！ 对于这个问题，我们可以先用欧拉公式， 对于这个问题，我们可以先用欧拉公式，通过给定的初值计 算出第一个点的值，然后在利用这两个（第一和第二个点） 算出第一个点的值，然后在利用这两个（第一和第二个点） 的值进行计算，直到计算出全部节点上的值。 的值进行计算，直到计算出全部节点上的值。 下面，我们用中点公式求解刚才的例子，计算的步长取 下面，我们用中点公式求解刚才的例子，计算的步长取0.01， ， 可以看到，计算的精度比较高。 可以看到，计算的精度比较高。 此时， 此时，计算公式为
常微分方程初值问题的的数值解法
本章讨论常微分方程初值问题的数值解法
2
考虑一阶常微分方程的初值问题
⎧ dy ⎪ = f ( x, y ) ⎨ dx ⎪ ⎩ y (a ) = y0
x ∈ [a, b]
只要 f (x, y) 在[a, b] × R1 上连续，且关于 y 满足 Lipschitz 条 件，即存在与 x, y 无关的常数 L 使对任意x∈[a, b] ，和y1， y2 ∈ R1 都有 | f ( x, y1) − f ( x, y2 ) | ≤ L| y1 − y2 | 在唯一解。 成立， 则上述问题存
⎧ ⎪ ⎨ ⎪ ⎩ y n +1 = yn + hf ( xn , yn ), h yn +1 = yn + [ f ( xn , yn ) + f ( xn +1 , y n +1 )] 2
改进的Euler方法：y0=1,
y1=y0+hf (x0, y0) =1.1, y1=1+01./2 ×[(1−2 ×0/1)+(1.1−2 ×0.1/1.1)] =1.095909, …… y11=…… y11=1.737869.
1 yn +1 = yn + h[ f ( xn , yn ) + f ( xn +1 , yn +1 )] 2
12
称之为梯形公式。这是一个隐式的计算公式，欲求的yn+1需 解一个方程。
3．截断误差
定义 在假设 yn = y(xn)，即第 n 步计算是精确的前提下，考 虑的截断误差 εn+1 = y(xn+1) − yn+1 称为局部截断误差
⎧ y n +1 = y n + k1 ⎨ ⎩k1 = hf ( xn ,y n )
常微分方程初值问题的数值解法
常微分方程初值问题的数值解法在实际应用中，对于某些微分方程，我们并不能直接给出其解析解，需要通过数值方法来求得其近似解，以便更好地理解和掌握现象的本质。
常微分方程初值问题（IVP）即为一种最常见的微分方程求解问题，其求解方法有多种，本文将对常微分方程初值问题的数值解法进行较为详细的介绍。
一、欧拉法欧拉法是最基本的一种数值解法，它采用泰勒级数展开并截断低阶项，从而获得一个差分方程近似求解。
具体来讲，设 t 为独立变量，y(t) 为函数 y 关于 t 的函数，方程为：$$y'(t) = f(t, y(t)), \qquad y(t_0) = y_0$$其中 f(t,y(t)) 为已知的函数，y(t_0) 为已知的初值。
将函数 y(t) 进行泰勒级数展开：$$y(t+h) = y(t) + hf(t, y(t)) + O(h^2)$$其中 h 表示步长，O(h^2) 表示其他高阶项。
为了使误差较小，一般取步长 h 尽可能小，于是我们可以用欧拉公式表示数值解：$$y_{n+1} = y_n + hf(t_n, y_n), \qquad y_0 = y(t_0)$$欧拉法的优点是容易理解和实现，但是由于截取低阶项且使用的单步法，所以误差较大，精度较低，在具体应用时需要慎重考虑。
二、龙格-库塔法龙格-库塔法（Runge-Kutta method）是一种多步法，比欧拉法更加精确。
龙格-库塔法的主要思想是使用不同的插值多项式来计算近似解，并且将时间步长分解，每次计算需要多次求解。
以下简要介绍二阶和四阶龙格-库塔法。
二阶龙格-库塔法将时间步长 h 分解成两步 h/2，得到近似解表达式：$$\begin{aligned} k_1 &= hf(t_n, y_n)\\ k_2 &= hf(t_n+h/2,y_n+k_1/2)\\ y_{n+1} &= y_n+k_2+O(h^3)\\ \end{aligned}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法，其需要计算的中间值较多，但是具有更高的精度。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
//输出计算结果
7
8.2 欧拉方法 显式欧拉法对应的程序
#include <stdio.h>
double f(double x,double y)
5
8.2 欧拉方法
一、显式欧拉法（续）
显式欧拉法的递推公式是斜率为f(xk-1,yk-1)，经过点(xk-1,yk-1)的直线方程。 上述显式欧拉法递推过程的几何含义，就是用曲线y(x)在点(x0,y0)处的切 线段代替y(x)在区间[x0,x1]内的曲线段；再把曲线段的终点(x1,y(x1))近似 为切线段的终点(x1,y1)，把曲线y(x)在点(x1,y1)处切线的斜率f(x1,y(x1))近 似为f(x1,y1) ，做曲线y(x)在点(x1,y1)处的切线，并在区间[x1,x2]内用近似 的切线段代替y(x)的曲线段；……。欧拉法用一系列折线近似地代替曲线 y(x)，因此它又称为欧拉折线法。
6
8.2 欧拉方法
显式欧拉法的算法
输入求解区间的边界a,b，步长h，起点x[0]=a处的纵坐标y[0]。
递推的次数n=(b-a)/h。
for(i=0;i<n;i++)
//循环1次完成1轮递推
x[i+1]=x[i]+h;
y[i+1]=y[i]+h*f(x[i],y[i]);
for(i=0;i<=n;i++) 输出x[i],y[i]。
在本章中假设讨论的一阶常微分方程的初值问题的解y(x)存在、唯一且 足够光滑，方程本身是稳定的，即精确解y(x)连续且依赖于初始值及右 端函数。
2
8.1 基础知识
二、数值解法
一阶常微分方程初值问题的数值解法的主要思想，是对区间[a,b]上的节点 a=x0＜x1＜……＜xn＜xn+1＜……≤b
建立y(xn)的近似值yn的某一递推格式，利用初值y0和已计算出的y1,y2,……, yk-1递推出yk，并且用这一方法反复递推，依次得到yk+1,yk+2,……,yn。这一 求解方法称为步进式求解，相邻2个节点的距离称为步长，记为hi=xi+1－xi。 为便于计算，常取成等距节点，称为定步长，这时把步长记为h。
由基于数值微分的离散化方法、基于数值积分的离散化方法、基于泰 勒展开的离散化方法都可以推导出显式欧拉法的递推公式。
用显式欧拉法求解一阶常微分方程初值问题的过程，就是以已知的(x0,y0) 作为起点，代入显式欧拉法的递推公式的右端，计算出y(x)在x1处的近似 值y1；再以(x1,y1)作为起点，用显式欧拉法的递推公式计算出y(x)在x2处 的近似值y2；……。
第8章 常微分方程初值问 题的数值解法
基础知识 欧拉方法 龙格—库塔方法
1
8.1 基础知识
一、问题的提出
很多实际问题都需要求解常微分方程。例如单摆问题。
常微分方程分为线性常微分方程和非线性常微分方程，又可以分为一阶 常微分方程和高阶常微分方程。通过变量的替换，可以把高阶常微分方 程转化为一阶常微分方程再求解。对于一阶常微分方程组，可以写成向 量形式的单个方程，求解方法与一阶常微分方程相似。因此本章只讨论 一阶常微分方程的初值问题：
dy f (x, y) , a x b dx y(x0 ) y0
目前在常微分方程理论中，只能求出某些特殊类型常微分方程的解析解， 对大部分常微分方程，用解析方法求出常微分方程的精确解非常困难， 甚至不存在解的解析表达式。为满足工程实践的需要，常常用数值解法 求常微分方程的近似解。
一阶常微分方程初值问题的数值解法有多种分类方法。 一种分类方法为： ⑴ 单步法：每一轮递推只用到前面一轮的递推结果，递推格式为： yk=yk-1＋hT(xk-1,yk-1) ⑵ 多步法：每一轮递推要用到前面多轮递推的结果，递推格式为： yk=yk-1＋hT(xk-r,yk-r,xk-r+1,yk-r+1,……, xk-1,yk-1)，其中r＞1。
定义：假设某一步递推的起点是精确的，即yi-1=y(xi-1)，那么这一步递推 的截断误差Ri=y(xi)－yi称为局部截断误差。 定义：若某算法的局部截断误差为O(hp+1)，则称此算法有p阶精度。
显式欧拉法具有一阶精度。
显式欧拉法在步长过大时误差较大；在步长较小时需要多步递推，可能 出现误差积累的现象。由于显式欧拉法的精度不高，因此在实际应用中 用的较少,……,yk-1，才能启动一 个r步的多步法。
3
8.1 基础知识
二、数值解法（续）
另一种分类方法为：
⑴ 显式方法：递推公式的右端都是已知量，可以直接计算出递推的结 果，递推格式为：yk=yk-1＋hT(xk-r,yk-r,xk-r+1,yk-r+1,……, xk-1,yk-1) ⑵ 隐式方法：递推公式左端的未知量也出现在公式的右端，递推格式 为： yk=yk-1＋hT(xk-r,yk-r,xk-r+1,yk-r+1,……, xk,yk) 隐式方法的递推公式其实是一个方程。解方程的运算量可能较大，为避免 解方程，常采用预测—校正系统。
4
8.2 欧拉方法
一、显式欧拉法
在一阶常微分方程初值问题的数值解法中，显式欧拉（Euler）法是最简 单的一种。显式欧拉法有明显的几何含义，缺点是精度不高。对于一阶 常微分方程的初值问题：
dy f (x, y) , a x b dx y(x0 ) y0
显式欧拉法的递推公式为：yk=yk-1＋hf(xk-1,yk-1)，k=1,2,3,……。 显式欧拉法每一轮递推只用到前面一轮递推的结果，因此它是单步法。
⑶ 预测—校正系统：每一轮递推包括预测和校正这2个步骤。先用显式方 法计算出yk，作为迭代的初值，这一过程称为预测；再把隐式方法的递推 公式作为迭代公式，把预测值yk代入迭代公式右端进行迭代，这一过程称 为校正。在校正时往往迭代1次或几次，校正值的精度就会有大幅提高。 一阶常微分方程初值问题的数值解法一般是对连续的初值问题进行离散化 处理，把微分方程转化为代数方程来求解。常用的离散化方法有： ⑴ 基 于数值微分的离散化方法，⑵ 基于数值积分的离散化方法，⑶ 基于泰勒 展开的离散化方法。