CFD2020-第10讲-常微分方程数值解法
第10章_常微分方程数值解
迭代法的收敛条件：
计
算 方
|
y y (k 1) i 1
(k ) i 1
|
h| 2
f
( xi1,
y(k) i 1
)
f
( xi1,
y (k 1) i 1
)
|
法 课
h 2
L|
y(k) i 1
y (k 1) i 1
|
件
其中L为Lipschitz常数。
当
0 h L 1
2
时迭代过程收敛。即迭代过程中误差没有扩大。
称为欧拉中点公式.容易看出,中点公式计算yi+1时,不仅需要yi的 值，还需要yi-1的值.
8
结束
凡是从已知(或已算出)的y0,y1,…,yi能直接从公式算出yi+1的公 式称显式，否则称隐式公式(8.4),(8.6)都是显式，(8.5)是隐式
在计算yi+1时，只需要yi的值，则称公式为单步法；若除yi之外 还需要以前的yi-1等多个值，则称多步法公式.(8.4)、(8.5)是单 计 步法，(8.6)是多步法，确切说是二步法.
0.2 xi yi
计 算
yi
1
yi
0.05
yi
2xi yi
y(0) i 1
2 xi 1 y (0)
i 1
方
y0 1
法 课
件 结果（略）
15
结束
§8.3 龙格-库塔(Runge-Kutta)法
欧拉方法是显式的一步法,使用方便,但精度较低.本节将构造出
高精度的显式一步法：龙格-库塔法,简称R-K法. 计 算 8.3.1 二阶R-K法
8.2.1 欧拉方法的导出
把区间［a,b］分为n个小区间，取步长h=(b-a)/ n ，节点 xi=x0+ih,i=0,1,2,…，n，其中x0=a，又设y (x)为上述问题的解.
常微分方程数值解法
用分段的折线逼近函数,此为 “折线法”而非“切线法”， 除第一个点是曲线上的切线，
其它都不是。
2、Euler方法的误差估计
1)局部截断误差。 在一步中产生的误差而非累积误差：
~
T x y y
n1
n1
n1
其中
~
y
是当
y
n
y(
x
)
n
（精确解！）时
n1
由Euler法求出的值，即y 无误差！ n
y y x , y x y 则得： h f
f ,
n1
n2
n
n
n1
n1
同样与Euler法结合，形成迭代算法，对n 0，1，2，
y y x , y 0 hf
n1
n
n
n
y y x , y x y
k 1
推出总体误差与步长的关系。
由微分方程解的存在唯一性,自然假定 （ f x，y）
充分光滑，或满足 Lipschitz条件：
f
x,ny源自xn
f
x
,
n
y
n
L
yxn
y n
第 n 步 的 总 体 截 断 误 差 记 为
en y
xn
y n
则 对 n 1 步：
e x y x y y y T y y ~ ~
用yn1, yn代替y(xn1), y(xn ), 对右端积分采用 取左端点的矩形公式
则有
xn1 xn
f
(x,
y)dx
hf
(xn ,
yn )
常微分方程的数值解法
常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。
然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。
在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。
这些方法可以给出解的近似表达式，通常称为近似解析方法。
还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。
利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题的精确解记为y (x )。
数值方法的基本思想是：在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。
这些h i 可以不相等，但一般取成相等的，这时na b h -=。
在这些节点上采用离散化方法，（通常用数值积分、微分。
泰勒展开等）将上述初值问题化成关于离散变量的相应问题。
把这个相应问题的解y n 作为y (x n )的近似值。
这样求得的y n 就是上述初值问题在节点x n 上的数值解。
一般说来，不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法 1.欧拉法1．对常微分方程初始问题(9.2))((9.1) ),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
常微分方程的数值解法
常微分方程的数值解法…………江南大学信计1203柯恒一、前言对于很多微分方程，我们很难求出解析解，这时我们需要采取数值手段求解。
在数值分析这门课中，老师讲到dy dt =f (t ,y )的数值解法，我们采用了欧拉格式，后退欧拉格式，梯形公式，改进欧拉格式及4阶龙格-库塔方法求解。
而老师没讲关于微分方程组和高阶微分方程的解法。
现在我来粗虐的说一下微分方程组及高阶微分方程的数值解，这里主要是借助matlab 中的相关函数进行计算并仿真出图像。
二、相关问题初值问题问题1，微分方程组问题描述：给定一个微分方程，并且告诉我们初始状态。
我们便可求出整个过程的解。
给定下面方程（L 表示省略号）：11122112112(,,,,)(,,,,)(,,,,)n n n n dy f x y y y dx dy f x y y y dx dy f x y y y dx ⎧=⎪⎪⎪=⎪⎨⎪⎪⎪=⎪⎩而且初始状态y ′i (x 0)已知记成y(0).问题解决：四阶龙格库塔方法：K i1=f i (x n ,y 1n ,y 2n ,……,y nn );K i2=f i (x n +h 2,y 1n +h 2∗k 11,y 2n +h 2∗K 21,……,y nn +h 2∗K n1) K i3=f i (x n +h ,y 1n +h ∗K 12,y 2n +h ∗K 22,……,y nn +h ∗K n2) K i4=f i (x n +h 2,y 1n +h 2∗K 13y 2n +h 2∗K 32,……,y nn +h 2∗K n3),y i,n+1=y i,n+h6∗(K i1+2∗K i2+2∗k i3+2∗K i4)通过龙格库塔方法，我们可以计算后面点的值，在matlab中调用ode45来实现四阶龙格库塔方法的调用。
应用举例：有一个同步地球卫星，现在加速到4km/s进行变轨，试分析变轨后卫星运动的轨迹。
已知地球的质量M=5.97e24,引力常量的值为6.672e-11问题建模：我们已经知道行星的运动是在一个平面上的。
数值分析常微分方程数值解
许多实际问题的数学模型是微分方程或微分方程的定解问题。
如物体运动、电路振荡、化学反映及生物群体的变化等。
常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。
若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。
因此研究一阶微分方程的初值问题⎪⎩⎪⎨⎧=≤≤=0)(),(y a y bx a y x f dxdy， (9-1) 的数值解法具有典型性。
常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。
用解析方法只能求出线性常系数等特殊类型的方程的解。
对非线性方程来说，解析方法一般是无能为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。
因此研究微分方程的数值解法是非常必要的。
只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。
由常微分方程的理论知，如果(9-1)中的),(y x f 满足条件（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续； （2）),(y x f 在上关于满足Lipschitz 条件，即存在常数，使得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称为由到1+n x 的步长。
今后如无特别说明，我们总假定步长为常量。
建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： (1) 用差商近似导数在问题(9-1)中，若用向前差商hx y x y n n )()(1-+代替)(n x y '，则得)1,,1,0( ))(,()()(1-=≈-+N n x y x f hx y x y n n n n n)(n x y 用其近似值代替，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=-这样，问题(9-1)的近似解可通过求解下述问题100(,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-⎧⎨=⎩(9-2)得到，按式(9-2)由初值经过步迭代，可逐次算出N y y y ,,21。
数值分析常微分方程数值解法
第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
常微分方程数值解法5262115页PPT文档
r 表示食饵独立生存时的增长率；
d 表示捕食者独立生存时的死亡率；
a 表示捕食者的存在对食饵增长的影响系数，反映捕
食者对食饵的捕获能力；
b 表示食饵的存在对捕食者增长的促进系数，反映食
饵对捕食者的喂养能力
150 100
令 y 1 y ,y 2 y ',y 3 y '', ,y n y ( n 1 )
可以将以上高阶微分方程化为如下一阶常微分方程组
y1 ' y2 y2 ' y3 yn ' an(x)y1
a1(x)yn f (x)
例：P120,1(a),Bessel方程
常微分方程的数值解
一般地，凡表示未知函数，未知函数的导 数与自变量之间的关系的方程叫做微分方 程．未知函数是一元函数的，叫常微分方 程;未知函数是多元函数的，叫做偏微分方 程．
如
y ' x y'x2y2 y''y'xy
Matlab实现 [t,x]=ode45(f,ts,x0,options,p1,p2,......)
50 0 0
30 20 10
0 0
10
20
50
30
20
10
0
30
0
10
8
6
4
2
100
0
50
100
150
50
100
高阶常微分方程的解法
高阶常微分方程
y ( n ) a 1 ( x ) y ( n 1 ) a ( n 1 ) ( x ) y ' a n ( x ) y f( x )
常微分方程的数值解法
第六章 常微分方程的数值解法 §6.0 引言§6.1 算法构造的主要途径§6.2 Runge-Kutta Method 算法§6.3 线性多步法§6.4 线性多步法的一般形式§6.5 算法的稳定性、收敛性§6.0 引 言1． 主要考虑如下的一阶常微分方程初值问题的求解:()()00,dy f x y dx y x y ⎧=⎪⎪⎨=⎪⎪⎩ 微分方程的解就是求一个函数y=y(x)，使得该函数满足微分方程并且符合初值条件。
2. 例如微分方程:xy '-2y=4x ；初始条件: y(1)=-3。
于是可得一阶常微分方程的初始问题24(1)3y y 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 称为微分方程的数值解。
如图所示:§6.1 算法构造的主要途径x 0 x 1 x 2 ...1 欧拉公式1.1 构造的思想：利用差商代替一阶导数，即010()()x x y x y x dy dx h =-≈，则 1000()()(,)y x y x f x y h -≈。
常微分方程数值解法
第七章 常微分方程数值解法常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。
另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。
因此，有必要探讨常微分方程初值问题的数值解法。
本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。
第一节 欧拉法求解常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy（1）的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<<<n x x x x 210 上的近似值 ,,,,,210n y y y y{}n y 称为问题的数值解，数值解所满足的离散方程统称为差分格式，1--=i i ix x h 称为步长，实用中常取定步长。
显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。
定理 设函数()y x f ,在区域+∞≤≤-∞≤≤y b x a D ,:上连续，且在区域D 内满足李普希兹(Lipschitz)条件，即存在正数L ，使得对于R 内任意两点()1,y x 与()2,y x ,恒有()()2121,,y y L y x f y x f -≤-则初值问题(1)的解()x y 存在并且唯一。
一、欧拉法(欧拉折线法)若将函数)xy 在点nx处的导数()n x y '用两点式代替， 即()hx y x y x y n n n )()(1-≈'+，再用n y 近似地代替()n x y ,则初值问题(1)变为⎩⎨⎧==++=+ ,2,1,0),()(001n x y y y x hf y y n n n n(2)(2)式就是著名的欧拉(Euler)公式。
以上方法称 为欧 拉法或欧拉折线法。
常微分方程数值解法
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 n是准确的，准确解y(xn 1 ) 和离散解y n 1之间的误差称为局部截断误差 记做： R n 1 y(xn 1 ) y n 1 下面我们分析Euler方法的局部截断误差，对于其他方法， 仅给出结论。
由于我们在局部截断误差中，假设y n是准确的，则y n y( xn ) 则格式：yn 1 yn hf ( xn , yn ) 可改写为： yn 1 y( xn ) hf ( xn , y( xn )) 另外我们可以得到Talor展开如下：
初值问题的Euler方法
类似地，我们可以得到：（自己回去推导） y ( xn 1 ) y ( xn ) f ( xn 1 , y ( xn 1 )) h 以yn 表示y ( xn )的近似值，则有 yn 1 yn hf ( xn 1 , yn 1 ) （n 0,1, 2,...） (2) 这就是隐式的`Euler公式或向后Euler方法，它与显式 的不同在于，它每算一步要解函数方程(2)才能得到 yn 1。
初值问题的Euler方法
如果取以上两式的算术 平均值的结果，则得 h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )] （n 0,1,2,... ） 2 称为梯形公式。 计算yn时常用以下迭代式：
(0) yn yn hf ( xn , yn ) 1 （3） h ( k 1) (k ) yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )] （k 0,1,2,... ） 2 ( k 1) (k ) 当 | yn y 1 n | ( k 1) 时, 取yn 1 yn 1
常微分方程数值解法ppt课件
若存在正的常数 L 使：
(Lipschitz)条件
|f( x ,y 1 ) f( x ,y 2 ) | L |y 1 y 2 | ( 1 .3 )
使 得 对 任 意 的 x [ a , b ] 及 y 1 ,y 2 都 成 立
则称 f (x,y) 对y 满足李普希兹条件，L 称为 Lipschitz常数.
节点 x i a i h i , 一 般 取 h i h ( ( b a ) / n ) 即 等 距
要计算出解函数 y(x) 在一系列节点
ax 0x 1x nb
处的近似值 y y(x ) i 完整版PPT课件i
16
yf(x,y) axb (1 .1 )
y(x 0) y0
(1 .2 )
对微分方程(1.1)两端从 xn到 xn1 进行积分
在大量的实际方程中出现的函数起码的连续性都 无法保证，更何况要求阶的导数
求解数值解
很多微分方程 根本求不到 问题的解析解！
重要手段。
完整版PPT课件
7
5.常微分方程数值解法的特点 常微分方程的数值解法常用来求近似解
根据提供的算法 通过计算机
数值解法得到的近似 解(含误差)是一个 离散的函数表.
便捷地实现
欧拉方法的导出把区间ab分为n个小区间步长为要计算出解函数yx在一系列节点iiyyx?iiixaihhhban?????一般取即等距节点处的近似值01naxxxb?????1iiihxx??nn等分001112yfxyaxbyxy????????对微分方程11两端从1nnxx?到进行积分11nnnnxxxxydxfxyxdx??????11nnxnnxyxyxfxyxdx?????右端积分用左矩形数值求积公式22baggxdxbagaba???????gxfxyx?令11nnnnxxnnfxyxnnyyfxyxh??????得x0x11nnnnnnyxyxhyxyhfxy??????1
常微分方程的数值解法全文
第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 龙格-库塔方法
龙格-库塔方法的原理
常用的数值解法包括欧拉方法、龙格-库塔方法、改进的欧拉方法、预估 校正方法和步进法等。
数值分析常微分方程数值解法
（Numerical Methods for Ordinary Differential Equations ）
中心差商近似导数
y( x1 )
y( x2 ) y( x0 ) 2h
y( x2 ) y( x0 ) 2h f ( x1 , y( x1 ))
x0
x1
x2
yi1 yi1 2h f ( xi , yi ) i 1, ... , n 1
假设 yi1 y( xi1 ), yi y( xi ) ,则可以导出 Ri y( xi1 ) yi1 O(h3 ) 即中点公式具有 2 阶精度。
只要f (x, y)在a,b R1上连续, 且关于 y 满足
Lipschitz 条件，即存在与 x, y无关的常数 L 使 | f (x, y1) f (x, y2 ) | L | y1 y2 |
对任意定义在 a,b上的 y1 x, y2 x都成立，
则上述IVP存在唯一解。
求函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值
yi y(xi ) (i 1, ... , n)
的方法称为微分方程的数值解法。 y1,L , yn 称为微分方程的数值解。
称节点间距 hi xi1 xi (i 0, ... , n 1)为步长， 通常采用等距节点，即取 hi = h (常数)。
三、初值问题的离散化方法
离散化方法的基本特点是依照某一递推公式， 按节点从左至右的顺序依次求出y(xi )的近似 yi 值(i 1, ... , n),取 y0。
常微分方程数值解算法
常微分方程数值解算法常微分方程是在物理、经济、生物、环境科学等领域中最基本的数学工具之一。
为了解决实际问题，需要求解这些方程的解。
但是，大部分常微分方程是无法求得解析解的，因此需要通过数值方法来求解。
在数值方法中，其基本思想是将微分方程化为一个逐步求解的问题。
通过离散化得到一个差分方程，然后通过数值方法求解这个差分方程。
本文将就常微分方程的数值解算法进行介绍和探讨。
1.欧拉方法欧拉方法是最基本的一种常微分方程数值解方法。
它的基本思想是将微分方程化为差分方程。
欧拉方法是一种一阶的显式方法。
通过计算当前点处的斜率即可进行逼近。
如下所示：y(t + h) = y(t) + hf(t, y(t))其中，h是步长。
f(t, y)是微分方程右边的函数。
欧拉方法的由来是其是以欧拉为名的。
这种方法的优点是简单明了，易于理解。
但是，其与真实解的误差随着步长增大而增大，误差不精，计算速度较慢等缺点也使其并非一个完美的数值解方法。
2.改进的欧拉方法改进的欧拉方法被认为是欧拉方法的一个进化版。
它是二阶数值方法，明显优于欧拉方法。
其基本思想是通过步长的平均值h/2来进行逼近。
y(t + h) = y(t) + h[ f(t, y(t)) + f(t + h, y(t) + hf(t, y(t))/2) ]其优点是能够更准确地逼近微分方程的解，只比欧拉方法多计算一些，但是其步长的误差随着步长增大而减小，并且计算速度比欧拉方法稍快。
因此，改进的欧拉方法是比欧拉方法更好的方法，效果相对较好。
3.龙格库塔方法龙格库塔方法是一种经典的数值解方法。
对于非刚性的方程可以得到较为精确的数值解。
其算法思路是利用多阶段迭代的方式，求解一些重要的插值点，并利用插值点的结果来逼近方程的解。
其公式如下：y(t + h) = y(t) + (h/6)*(k1 + 2k2 + 2k3 + k4)其中，k1 = f(t, y(t))k2 = f(t + h/2, y(t) + h/2k1)k3 = f(t + h/2, y(t) + h/2k2)k4 = f(t + h, y(t) + hk3)其优点是更精确，计算速度更快。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
t x
y
z
x
y
z
FDM
空间离散
U Q(U)
FVM
时间推进 t
DG
推荐方法： Runge-Kutta法
2阶
U (1) U n tQ(U n )
U n1 1/ 2U n 1/ 2U (1) 1/ 2tQ(U (1) )
3 阶 （TVD型）
U (1) U n tQ(U n ) U (2) 3/ 4U n 1/ 4U (1) 1/ 4tQ(U ) (1) U n1 1/ 3U n 2 / 3U (2) 2 / 3tQ(U ) (2)
仅迭代一次(n=1),将
y(0) m 1
写入第二式
ym1
ym
h 2
f (xm , ym ) f (xm1, ym hf (xm , ym ))
改进的Euler方法
Copyright by Li Xinliang
7
P阶Taylor级数法
设初值问题的解具有p+1阶连续导数，利用Taylor公式
Euler公式
K1 f (xm , ym )
ym1 ym hf (xm , ym )
N 2
ym1 ym h b1K1 b2K2
K1 f (xm , ym )
K2 f (xm c2h, ym ha21K1)
• 将 K2 在(xm , ym ) 处展开
K2
f
h(c2
fx
a21K1 f y )
y'''(xm ) O(h4 )
y' f
利 用
y'' f x ff y
y''' fxx 2 ffxy f 2 f yy ( fx ffy ) f y
二阶：O(h3 )
b1 b2 1 b2c2 1/ 2
c2 a21
无穷多组解
b2a21 1/ 2
Copyright by Li Xinliang
ym1 ym h(xm , ym , h)
d dx
f
(x,
y(x))
fx
f fy
d2 dx2
f (x, y(x))
f xx 2 ff xy
f 2 f yy
fy( fx
ff y )
……
(*)式导数的计算比较复杂 离散情形（非解析形式）
（很少使用）
Copyright by Li Xinliang
向前Euler方法
ym1 ym hf (xm1, ym1)
向后Euler方法
ym1
ym
h 2
f (xm , ym ) f (xm1, ym1)
梯形法
迭求代解
y(0) m1
y(n) m1
ym ym
hf (xm , ym )
h 2
f (xm , ym )
f (xm1,
y ) (n1) m1
K1 f (xm , ym )
K2 f (xm c2h, ym ha21K1)
f
h(c2
fx
a21K1 f y )
h2 2
(c22
f xx
2c2a21K1 f xy
a221K12
fyy)
h3 6
(c23 f xxx 3c22a21K1 f xxy 3c2a221K12
f
xyy
y(xm1)
p j0
hj j!
y
(
j)
(
xm
)
(
h p
p1
1)!
y(
p1)
(m
)
y' f (x, y)
(x, y(x),h)
p h j1 j1 j!
d j1 dx j1
f (x, y(x))
(*)
y ( xm 1 )
y ( xm
)
h(xm ,
y(xm ),
h)
(
h p1 p 1)!
y( p1) (m )
f yyy)
m
O(h4 )
ym1 ym b1hK1 b2hK2 b3hK3 O(h4 )
ym b1hf b3h[ f
b2h[ h(c3 fx
f
h(c2 fx (a31K1 a32
a21K1 K2) fy
f )
y)
h2
2
h2 2
(c22
(c32 f xx
f xx 2c2a21K1 f xy 2c3 (a31K1 a32K2
• ym1 ym h(b1 b2 ) f h2b2 (c2 fx a21 ffy )
h3 2
b2 (c22
f xx
2c2a21
ff
xy
a221
f
2
f yy)
O(h4 )
• 将 y(xm1) 在 xm 处展开
y ( xm1 )
y(xm ) hy'(xm )
h2 2
y''(xm )
h3 3!
Rm1 y(xm1) y(xm ) h(xm , y(xm ), h)
从准确值y( xm )计算1步后得 到的值与准确值y(xm1) 之差
整体截断误差:
m1 y(xm1) ym1
P阶方法:
从初值出发，计算m+1步后 得到的值与准确值之差
局部截断误差 O(h p1) 显式方法:
右端项不包含 ym1 隐式方法:
xm1 xm
f (x, y(x))dx
• 利用Taylor公式
y ( xm1 )
y(xm ) hy'(xm )
h2 2
y''(m )
Copyright by Li Xinliang
5
单步法
计算ym1 的值只用到前一步的值 ym
ym1 ym h(xm , ym , h)
增量函数
局部截断误差:
h3 2
b2 (c22
f xx
2c2a21
ff
xy
a221
f
2
f yy)
O(h4 )
Copyright by Li Xinliang
10
Runge-Kutta方法
N阶
N
ym1 ym h bi Ki i 1
N 2 ym1 ym h b1K1 b2K2
K1 f (xm , ym ) K2 f (xm c2h, ym ha21K1)
4阶
U (1) U n 1/ 2tQ(U n ) U (2) U n 1/ 2tQ(U ) (1) U (3) U n tQ(U (2) ) U n1 1/ 3(U n U (1) 2U (2) U ) (3) 1/ 6tQ(U ) (3)
更高阶 ……
Copyright by Li Xinliang
右端项含有 ym1
• 一阶向后差商代替一阶导数
y(xm1) h
y(xm )
y' ( xm1 )
• 两边从xm 到xm1 积分,积分采用 右矩形公式
y(xm1) y(xm )
xm1 xm
f (x, y(x))dx
Copyright by Li Xinliang
6
改进的Euler方法
ym1 ym hf (xm , ym )
y
y' f (x, y)
y2
y1
y(x0 ) y0
初值条件
y0
因为f(x,y(x))恰好是待求精确解y(x)的斜率y’(x)
y(x)是在x-y平面上经过点(x0,y0)的曲线,该曲线
在点(x0,y0)处的切线的斜率为f(x0,y0).
y y(x)
x0
x1 x2
x
Euler方法
y y0 f (x0 , y0 )( x x0 )
x x1 x h
y1 y0 hf (x0 , y0 ) y2 y1 hf (x1, y1)
ym1 ym hf (xm , ym )
• 一阶向前差商代替一阶导数
y(xm1) h
y(xm )
y'(xm )
• 两边从 xm到xm1 积分,积分采用
左矩形公式
y(xm1) y(xm )
11
Runge-Kutta方法
N阶
N 2
ym1 ym h b1K1 b2K2
K1 f (xm , ym ) K2 f (xm c2h, ym ha21K1)
2阶格式满足的条件：
N
ym1 ym h bi Ki i 1
b1 b2 1 b2c2 1/ 2 b2a21 1/ 2
y Ce f (x)dx
y(x) c'(x) g(x)e f (x)dx
• 常微分方程的数值计算
P(x, y) Q(x, y) y' 0
U P, x
U Q y
d U (x, y(x)) 0 dx
U (x, y(x)) C
计算机的发展、绝大多数常微分方程难以求解析解
U f1(U) f2 (U) f3 (U) g1(U) g2 (U) g3 (U)
f (x, y) 1
(y x)
i 1
ci aij , i 2,, N j 1
新的ym1与原来的 x 在计 算过程中有相同的增量
N
bi 1
i 1
Copyright by Li Xinliang
13
N
Runge-Kutta方法
ym1 ym h bi Ki
i 1
N 3 将 K2 K3 在 (xm, ym ) 处展开
K1 f (xm , ym )
i 1
Ki f (xm cih, ym h aij K j ) j 1