最新9常微分方程数值解59382
研常微分方程的数值解法.pptx
y0
a x b 一阶方程初值问题
y y(a)
f
( x, y0
y, y) a x , y(a) y0
b
高阶方程初值问题
y1
y2
f1( x, y1 , y2 ) f2 ( x, y1 , y2 )
y1( x0 ) y1 方程组初值问题 y2( x0 ) y2
.微分方程"解析解"存在的条件 P186 定理1
f ( xn1 , yn1 )
yn1
2 xn1 yn1
yn hK1
2 xn1 yn hK1
K2
h yn1 yn 2 (K1 K2 )
第10页/共30页
-10-
三 局部截断误差和方法的阶数 局部截断误差 将方程精确解y( x)代入数值求解公式左右两端， 左右两端之差T =y( xn1 ) yn1为方法的局部截断误差。 方法的精度 若T O(hp1 )，则称此方法具有p阶精度或称方法是p阶的。 Euler 公式的局部截断误差与精度
改进的Euler公式
-9-
例2
求解初值问题
y y 2x y
0 x1 （h 0.1）
y(0) 1
解 f (x, y) y 2x
y
Euler公式
yn1
yn
hf
( xn ,
yn )
yn
h( yn
2 xn yn
)
改进Euler公式
f
( xn ,
yn )
yn
2 xn yn
K1
yn1 yn hf ( xn , yn ) yn hK1
xn )
y( ) ( x 2!
xn )2
令x xn1 :
数值分析第九章常微分方程数值解法
松弛法
通过迭代更新函数值并逐步放松约束 条件来逼近解，适用于刚性和非刚性 问题。
利用线性组合迭代函数值来逼近解， 具有更高的收敛速度和稳定性。
03
数值解法的稳定性分析
数值解法的稳定性定义
数值解法的稳定性是指当微分方程的初值有微小的扰动时， 其数值解的近似值的变化情况。如果数值解在微小扰动下变 化较小，则称该数值方法是稳定的。
更高的精度和稳定性。
数值逼近法
泰勒级数法
将微分方程的解展开为泰勒级数，通过截断级数来逼 近解。
多项式逼近法
利用多项式来逼近微分方程的解，通过选取合适的基 函数和系数来提高逼近精度。
样条插值法
利用样条函数来逼近微分方程的解，具有更好的光滑 性和连续性。
迭代法
雅可比迭代法
通过迭代更新函数值来逼近微分方程 的解，具有简单易行的优点。
初值和边界条件的处理
根据实际问题，合理设定初值和边界 条件，以获得更准确的数值解。
收敛性和误差分析
对数值解进行收敛性和误差分析，评 估解的精度和稳定性。
数值解法的应用案例分析
人口增长模型
通过数值解法求解人口增长模型，预测未来人口数量，为政策制 定提供依据。
化学反应动力学
利用数值解法研究化学反应的动力学过程，模拟反应过程和结果。
数值分析第九章常微分方 程数值解法
• 引言 • 常微分方程数值解法的基本思想 • 数值解法的稳定性分析 • 数值解法的收敛性和误差分析 • 数值解法的实现和应用案例
01
引言
常微分方程的应用背景
自然科学
描述物理、化学、生物等自然 现象的变化规律。
工程领域
控制系统设计、航天器轨道计 算等。
数值分析常微分方程数值解法
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
数值分析课件第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
常微分方程数值解法教材
常微分方程数值解法科学技术中的许多问题，在数学中往往归结为微分方程的求解问题。
例如天文学中研究星体运动，空间技术中研究物体飞行等，都需要求解常微分方程初值问题。
最简单的常微分方程初值问题是座=f (x, y), a £ x £ b,dx ' " (*)y(x o) =y°.定理如果f (x, y)在区域D ={( x,y) a《x〈b, y<°°}内连续，且关于y满足Lipschitz条件，那么初值问题 (*)的解存在且惟一。
因为数值解是求微分方程解y(x)的近似值，所以总假定微分方程的解存在惟一，即初值问题(*)中的f(x, y)满足定理的条件。
除特殊情形外，初值问题(*) 一般求不出解析解，即使有的能求出解析解，其函数表示式也比较复杂，计算量比较大，况且实际问题往往只要求在某一时刻解的函数值，因此，有必要研究初值问题(*)的数值解法。
所谓初值问题(*)的数值解法，就是寻求解y(x)在区间［a,b］上的一系列点x〔：：：x^：： x3 ：：：..• ::：x n:：:■■-上的近似值y i,y2,…，y n,….记h j = x —x口(i =1,2, |||)表示相邻两个节点的间距，称为步长。
求微分方程数值解的主要问题：(1)如何将微分方程也=f (x, y)离散化，并建立求其数值解的递推公式；dx(2)递推公式的局部截断误差、数值数y n与精确解y(x n)的误差估计；(3)递推公式的稳定性与收敛性.、Euler方法Euler法是最早的解决一阶常微分方程初值问题的一种数值方法，虽然它不够精确，很少被采用，但是它在某种程度上反映了数值方法的基本思想和特征。
考虑初值问题% f(x, y), (1.1)dxyMi。
，(1.2) 为了求得它在等距离散点x1 <x2<x n上的数值解，首先将(1.1)离散化。
设h =X i —x — (i =1,2,|||)，将式(1.1)离散化的办法有 Taylor 展开法、数值微分法及数值积分法。
常微分方程数值解法
i.常微分方程初值问题数值解法i.1 常微分方程差分法考虑常微分方程初值问题：求函数()u t 满足(,), 0du f t u t T dt=<≤ （i.1a ） 0(0)u u = (i.1b)其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。
我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-∀∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。
通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。
本章讨论常微分方程最常用的近似数值解法--差分方法。
先来讨论最简单的Euler 法。
为此，首先将求解区域[0,]T 离散化为若干个离散点：0110N N t t t t T -=<<<<= （i.3） 其中n t hn =，0h >称为步长。
在微积分课程中我们熟知，微商（即导数）是差商的极限。
反过来，差商就是微商的近似。
在0t t =处，在（i.1a ）中用向前差商10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到1000(,)u u hf t u -=一般地，我们有1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
下面我们用数值积分法重新导出 Euler 法以及其它几种方法。
为此，在区间1[,]n n t t +上积分常微分方程（i.1a ），得11()()(,())n n t n n t u t u t f t u t dt ++=+⎰ （i.5）用各种数值积分公式计算（i.5）中的积分，便导致各种不同的差分法。
常微分方程数值解-PPT精品文档
称为局部截断误 差。显然，这个 y ( x ) y ( x ) h 误差在逐步计算 n 1 n y ' ( x ) y ' ' ( ) n n 过程中会传播， h 2 积累。因此还要 y ( x ) y ( x ) h n 1 n f ( x , y ( x )) y ' ' ( ) 估计这种积累 n n n h 2
对于一个常微分方程：
9.1 Euler方法
dy y ' f( x ,y ), x [ a , b ] dx 通常会有无穷个解。如：
dy cos( x ) y sin( x ) a , a R dx 因此，我们要加入一个限定条件。通常会在端点出给出， 如下面的初值问题： dy f (x , y) , x [a ,b ] dx )y 0 y(a 为了使解存在唯一，一般，要加限制条件在f上，要求f对y 满足Lipschitz条件：
求 y ( x ) 在 x i 上的近似值
y i 。 { y i } 称为分割 I
上的格点函数
我们的目的，就是求这个格点函数
② 由微分方程出发，建立求格点函数的差分方程。这个方程应该满足： A、解存在唯一；B、稳定，收敛；C、相容 ③ 解差分方程，求出格点函数
数值方法，主要研究步骤②，即如何建立差分方程，并研究 差分方程的性质。
x0
x1
y i 1 y i h f ( x i 1 , yi 1 ) ( i 0, ... , n 1)
由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。 一般先用显式计算一个初值，再迭代求解。
常微分方程数值解法
而
y ( xn 1 ) y ( xn ) hy( xn ) 1 2 1 1 h y( xn ) h3 y ( xn ) y (4) ( xn ) O(h 5 ) 2! 3! 4!
要使y(xn+1)–yn+1= O(h5) ，比较h的同次幂系数，必 有
h : a 1 1 h : a b c d 1 2 h : a 2b 2d 1 h3 : a 3b 3d 1 h 4 : a 4b 4d 1
yn 1 yn hf ( xn , yn ) h yn 1 yn 2 [ f ( xn , yn ) f ( xn 1 , yn 1 )] (n 0,1, 2,)
局部截断误差为O(h3)，是二阶方法.
标准四阶Runge-Kutta公式
1 yn 1 yn ( K1 2 K 2 2 K 3 K 4 ) K hf ( x ,6y ) n n 1 h K1 ) K 2 hf ( xn , yn 2 2 K hf ( x h , y K 2 ) n n 3 2 2 K 4 hf ( xn h, yn K 3 )
用Euler预报—校正格式求解初值问题
dy y y 2 sin x 0 dx y (1.2)及y(1.4)的近似值， 小数点后至少保留5位. 解 设f(x,y)=-y-y2sinx , x0=1,y0=1, xn=x0+nh=1+0.2n, Euler预报—校正格式为
h yn 1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 2 0.2 yn 1 yn [8 3 yn 8 3 yn 1 ] 2
数值分析--第9章常微分方程数值解
数值分析--第9章常微分⽅程数值解数值分析--第9章常微分⽅程数值解第九章常微分⽅程数值解法许多实际问题的数学模型是微分⽅程或微分⽅程的定解问题。
如物体运动、电路振荡、化学反映及⽣物群体的变化等。
常微分⽅程可分为线性、⾮线性、⾼阶⽅程与⽅程组等类；线性⽅程包含于⾮线性类中，⾼阶⽅程可化为⼀阶⽅程组。
若⽅程组中的所有未知量视作⼀个向量，则⽅程组可写成向量形式的单个⽅程。
因此研究⼀阶微分⽅程的初值问题=≤≤=0)(),(y a y b x a y x f dx dy ， (9-1)的数值解法具有典型性。
常微分⽅程的解能⽤初等函数、特殊函数或它们的级数与积分表达的很少。
⽤解析⽅法只能求出线性常系数等特殊类型的⽅程的解。
对⾮线性⽅程来说，解析⽅法⼀般是⽆能为⼒的，即使某些解具有解析表达式，这个表达式也可能⾮常复杂⽽不便计算。
因此研究微分⽅程的数值解法是⾮常必要的。
只有保证问题(9-1)的解存在唯⼀的前提下，研究其数值解法或者说寻求其数值解才有意义。
由常微分⽅程的理论知，如果(9-1)中的),(y x f 满⾜条件（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续；（2）),(y x f 在D 上关于y 满⾜Lipschitz 条件，即存在常数L ，使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟⼀的连续解)(x y y =。
在下⾯的讨论中，我们总假定⽅程满⾜以上两个条件。
所谓数值解法，就是求问题(9-1)的解)(x y y =在若⼲点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的⽅法。
),,2,1(N n y n =称为问题(9-1)的数值解，n n x x h -=+1称为由n x 到1+n x 的步长。
今后如⽆特别说明，我们总假定步长为常量。
常微分方程的数值解法全文
第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 <<.。
常微分方程的数值解
欧拉方法的实现
确定步长和初始值
根据问题的性质和精度要求，选择合适的步长 和初始值。
迭代计算
根据欧拉方法的公式，迭代计算下一个点的值。
终止条件
当达到预设的迭代次数或误差范围时，停止迭代。
常微分方程的应用
总结词
常微分方程在自然科学、工程技术和社会科学等领域有广泛应用。
详细描述
在物理学中，常微分方程可以描述物体的运动规律、电磁波的传播等；在化学中，可以描述化学反应 的动力学过程；在社会科学中，可以用于研究人口增长、经济趋势等。此外，常微分方程还在控制工 程、航空航天等领域有广泛应用。
确定步长和初始值
在应用龙格-库塔方法之前，需要 选择合适的步长和初始值。步长 决定了迭代的精度，而初始值则 用于启动迭代过程。
编写迭代公式
根据选择的步长和初始值，编写 龙格-库塔方法的迭代公式。该公 式将使用已知的函数值和导数值 来计算下一步的函数值。
迭代求解
按照迭代公式进行迭代计算，直 到达到所需的精度或达到预设的 最大迭代次数。
欧拉方法的误差分析
截断误差
由于欧拉方法只使用了微分方程的一次项， 因此存在截断误差。
全局误差
全局误差是实际解与近似解之间的最大偏差。
局部误差
由于每一步都使用了上一步的结果，因此存 在局部误差。
稳定性
欧拉方法是稳定的，但步长和初始值的选择 会影响其稳定性和精度。
04 龙格-库塔方法
龙格-库塔方法的原理
常用的数值解法包括欧拉方法、龙格-库塔方法、改进的欧拉方法、预估 校正方法和步进法等。
第九章常微分方程数值解 (3)PPT课件
3．公式的截断误差 二元泰勒公式：
设 z=f(x,y) 在点 (x0, 的y0某) 一邻域内连续且直到有n+1阶 连续偏导数，(x0h,y0为此h)邻域内任一点，则有：
f
(x0
h, y0
k)
f
(x0
,
y0
)
h
x
k
y
f
(
x0
,
y0
)
2
21! h
x
k
y
f (x0 , y0 ) ...
准确解
1 1.095445 1.183216 1.264911 1.341641 1.414214 1.483240 1.549193 1.612452 1.673320 1.732051 1.788854 1.843909 1.897367 1.949359 2.000000
定义 在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考 虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差
定义 若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。
Ri 的主项
欧拉法的局部截断误差：
R i y ( x i 1 ) y i 1 [ y ( x i ) h y ( x i ) h 2 2 y ( x i ) O ( h 3 ) [ y i ] h ( x i , y i ) f]
第九章 常微分方程数值解
考虑一阶常微分方程的初值问题
dy f(x, y) x[a,b] dx y(a) y0
只要 f (x, y) 在[a, b] R1 上连续，且关于 y 满足 Lipschitz 条 件，即存在与 x, y 无关的常数 L 使 |f(x ,y 1 ) f(x ,y 2 )| L |y 1 y 2 | 对任意定义在 [a, b] 上的 y1(x) 和 y2(x) 都成立，则上述问题存 在唯一解。
第九章 常微分方程数值解法
1.10000000000000
1.19181818181818 1.27743783371472 1.35821259956029 1.43513291865780
1.09073753683521
1.17407576129348 1.25124850679651 1.32309349775209 1.39017807462719
9.0 引言
2.使用数值积分的方法 对于微分方程 y’=f(x,y)两边求x0到x的定积分
x dy x x0 dx dx x0 f ( x , y( x ))dx
x y( x ) y( x0 ) f ( t , y( t ))dt x0
即 或写为
x y ( x) y0 f (t , y (t ))dt x0
x0 0, N 10, b 1, y0 1
由Euler 公式,有
2 xn ) yn1 yn hf ( xn , yn ) yn h( yn yn
n 1, 2, ,10
2019/2/26
9.1 Euler方法
9.1.1 Euler方法
得
2 x0 y1 y0 h( y0 ) y0
y( xn1 ) y( xn )
xn1 xn
f (t , y(t ))dt
xn 1 xn f ( xn , y( xn )) f ( xn1 , y( xn1 )) 2
故其数值公式为
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )] 2
数值解
x y
x0 x1 y0 y1
x2 … xn … y2 … yn …
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
0
0.44
可以看出误差随着计算在积累。
➢ 欧拉公式的改进：
隐式欧拉法 /* implicit Euler method */
向后差商近似导数
y(x1)y(x1) hy(x0) x0
y( x1 ) y0 h f ( x1, y( x1 ))
§1 Euler’s Method x1
yi1 yi h f ( xi1 , yi1 ) (i 0, ... , n 1)
§1 Euler’s Method
梯形公式 /* trapezoid formula */ — 显/隐式两种算法的平均
y i 1 y i h 2 [f( x i,y i) f( x i 1 ,y i 1 )] ( i 0 ,.,.n .1 )
注：的确有局部截断误差 R iy(x i 1)yi 1 O (h 3)， 即梯形公式需具要有22个阶初精值度y，0和比y欧1来拉启方动法递有推了进步。 但注意过到程该，公这式样是的隐算式法公称式为，双计步算法时/* 不do得ub不le-用ste到p 迭代法m。ethod */，而前面的三种算法都是单步法 /* single-step method */。
中点欧拉公式 /* midpoint formula */
中心差商近似导数
y(x1)y(x2)2hy(x0)
y (x 2 ) y (x 0 ) 2 h f(x 1 ,y (x 1 ))
x0
x1
x2
y i 1 y i 1 2 h f ( x i,y i)i 1பைடு நூலகம்,.,. n . 1
假设 yi 1y(x i 1),yiy(x i)，则可以导出 R iy(xi 1)yi 1O (h 3) 即中点公式具有 2 阶精度。
由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。 一般先用显式计算一个初值，再迭代求解。
隐式欧拉法的局部截断误差：
R i y(xi1)yi1 h 22 y(xi)O (h3)
即隐式欧拉公式具有 1 阶精度。
例2：dy y2x, y(0)1
dx
y
解：梯形公式为：yk1ykh 2 yk2y x kk yk12y x kk 11
解yk+1出来比较困难，遇到的是一个二次方程，
y 2 k 1 (1 h 2 )y k 1 ((1 h 2 )y k h y x k k) h x k 1 0
改进欧拉法 /* modified Euler’s method */
§1 Euler’s Method
Step 1: 先用显式欧拉公式作预测，算出 yi1 yi h f ( xi , yi )
Step 2: 再将 yi1 代入隐式梯形公式的右边作校正，得到
yi 1
yi
h 2
[
f
(
x
i
,
yi
y5= y4+hf(x4, y4)=-0.3333 y6= y5+hf(x5, y5)=0
xk
y(xk)
yk
ek
1.2
- 0.96
-1
0.04
1.4
-0.84 -0.9333 0.0933
1.6
-0.64
-0.8
0.16
1.8
-0.36
-0.6
0.24
2.0
0
-0.3333 0.3333
2.2
0.44
xk
yk
y(xk)
ek
0.1
1.0959691
1.0954451
0.000464
0.2
1.1840966
1.1832160
0.000881
0.3
1.2662014
1.2649111
0.001290
0.4
1.3433602
1.3416408
0.001719
0.5
1.4164019
1.4142136
0.002188
)
f ( xi1,
yi1 )]
y i 1 y i h 2 f( x i,y i) fx i 1 ,y i h f( x i,y i) ( i 0 ,.,n . 1 .)
注：此法亦称为预测-校正法 /* predictor-corrector method */。 可以证明该算法具有 2 阶精度，同时可以看到它是个单 步递推格式，比隐式公式的迭代求解过程简单。后面将 看到，它的稳定性高于显式欧拉法。
方法 显式欧拉 隐式欧拉 梯形公式
中点公式
简单 稳定性最好 精度提高
精度提高, 显式
§1 Euler’s Method
精度低 精度低, 计算量大 计算量大
多一个初值, 可能影响精度
withCaaollnft’ththeyiDetoaodpudmioOyvWsgsoaaaKsiudvkenib,levetltlal,ahmeeigcntt?iena’etssalkalgymefeoste?rwmgirtuehleaoduyt…any possible.
9常微分方程数值解59382
例1：
d d
y x
2y x
2
y 1 1
步长h=0.2
解：h=0.2 , xi=1+ih， 精确解为：y=x2-2x
y1= y0+hf(x0, y0)=-1
y2= y1+hf(x1, y1)=-0.9333
y3= y2+hf(x2, y2)=-0.8 y4= y3+hf(x3, y3)=-0.6
0.6
1.4759556
1.4832397
0.002716
0.7
1.5525141
1.5491933
0.003321
0.8
1.6164748
1.6124515
0.004023
0.9
1.6781664
1.6733201
0.004746
1.0
1.7378674
1.7320508
0.005817
§3 龙格 - 库塔法 /* Runge-Kutta Method */
可即以：预测y 校p 正y k方 法h ( y来k 求2：yxky kk )1 y k h 2f(x k,y k) fx k 1 ,y k hf(x k,y k)
y
c
yk
h( y p
2 xk1 ) yp
yk1
1 2
(
y
p
yc)
y k 1 y k h 2 f(x k,y k) fx k 1 ,y k h f(x k,y k)
建立高精度的单步递推格式。
单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜 率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所 能达到的最高精度为2阶。