最新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
- 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阶。