数值分析常微分方程数值解
数值分析第九章常微分方程数值解法
松弛法
通过迭代更新函数值并逐步放松约束 条件来逼近解，适用于刚性和非刚性 问题。
利用线性组合迭代函数值来逼近解， 具有更高的收敛速度和稳定性。
03
数值解法的稳定性分析
数值解法的稳定性定义
数值解法的稳定性是指当微分方程的初值有微小的扰动时， 其数值解的近似值的变化情况。如果数值解在微小扰动下变 化较小，则称该数值方法是稳定的。
更高的精度和稳定性。
数值逼近法
泰勒级数法
将微分方程的解展开为泰勒级数，通过截断级数来逼 近解。
多项式逼近法
利用多项式来逼近微分方程的解，通过选取合适的基 函数和系数来提高逼近精度。
样条插值法
利用样条函数来逼近微分方程的解，具有更好的光滑 性和连续性。
迭代法
雅可比迭代法
通过迭代更新函数值来逼近微分方程 的解，具有简单易行的优点。
初值和边界条件的处理
根据实际问题，合理设定初值和边界 条件，以获得更准确的数值解。
收敛性和误差分析
对数值解进行收敛性和误差分析，评 估解的精度和稳定性。
数值解法的应用案例分析
人口增长模型
通过数值解法求解人口增长模型，预测未来人口数量，为政策制 定提供依据。
化学反应动力学
利用数值解法研究化学反应的动力学过程，模拟反应过程和结果。
数值分析第九章常微分方 程数值解法
• 引言 • 常微分方程数值解法的基本思想 • 数值解法的稳定性分析 • 数值解法的收敛性和误差分析 • 数值解法的实现和应用案例
01
引言
常微分方程的应用背景
自然科学
描述物理、化学、生物等自然 现象的变化规律。
工程领域
控制系统设计、航天器轨道计 算等。
数值分析9-2
n 0 1 2 3
xn 0 0.02 0.04 0.06
yn 1.0000 0.9820 0.9650 0.9489
y(xn) 1.0000 0.9825 0.9660 0.9503
n = y(xn) - yn
0 0.0005 0.0005 0.0014
1.483240 1.549193 1.612452 1.673320 1.732051
10 1
9.2.4． 单步法的局部截断误差与阶
显式单步法一般形式为 yn 1 yn h ( xn , yn , h) 而隐式单步法一般形式为 yn 1 yn h ( xn , yn , xn 1 , yn 1 , h) 函数与f ( x, y )有关，称为增量函数。
f ( x, y ) , x [ a , b ] dx y ( a ) y0 为了使解存在唯一，一般，要加限制条件在f上，要求f对y 满足Lipschitz条件：
f ( x, y1 ) f ( x, y2 ) L y1 y2
常微分方程的解是一个函数，但是，计算机没有办法对函 数进行运算。因此，常微分方程的数值解并不是求函数的近似， 而是求解函数在某些节点的近似值。 要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值 yi y( xi ) ( i 1, ... , n)
为了考察数值方法提供的数值解，是否有实用价值， 需要知道如下几个结论： ① 步长充分小时，所得到的数值解能否逼近问题的 真解；即收敛性问题 ② 误差估计
③产生得舍入误差，在以后得各步计算中，是否 会无限制扩大；稳定性问题
数值分析 第9章 常微分方程初值问题数值解法
9 .2 .2 梯形方法/* trapezoid formula */— 显、隐式两种算法的平均 为得到比欧拉法精度高的计算公式, 在等式( 2.4) 右端积分 中若用梯形求积公式近似, 并用yn 代替y ( xn ) , yn+1 代替y ( xn+1 ) ,则得
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )], 2
yn 1 yn f ( xn , yn ), xn 1 xn
即 yn+1 = yn + hf ( xn , yn ) . ( 2 .1)
这就是著名的欧拉( Euler ) 公式.
• 若初值y0 已知, 则依公式( 2.1)可逐步算出
• y1 = y0 + hf ( x0 , y0 ) ,
为了分析迭代过程的收敛性, 将( 2. 7) 式与(2. 8 )式相减, 得
h ( k 1) (k ) yn 1 yn [ f ( x , y ) f ( x , y 1 n 1 n 1 n 1 n 1 )] 2
于是有
| yn 1 y
( k 1) n 1
hL (k ) | | yn 1 yn 1 |, 2
| f ( x, y1 ) f ( x, y2 ) | L | y1 y2 |, y1, y2 R,
定理1 设f在区域Ｄ＝{(x,y)|a≤x ≤b,y∈R}上连续， 关于y满足利普希茨条件，则对任意x0 ∈[a,b], y0 ∈R，常微分方程初值问题(1.1)式和(1.２)式当x ∈[a,b]时存在唯一的连续可微解y(x). 定理２ 设f在区域Ｄ上连续，且关于y满足利 普希茨条件，设初值问题
1 2 1 2 dy x ydy xdx y x c 2 2 dx y y (0) 2 y2 x2 4
第9章 常微分方程初值问题数值解法
数值分析
第9章 常微分方程初值问题数值解法
《常微分方程》中介绍的微分方程主要有:
(1)变量可分离的方程 (2)一阶线性微分方程(贝努利方程) (3)可降阶的一类高阶方程 (4)二阶常系数齐次微分方程 (5)二阶常系数非齐次微分方程 (6)全微分方程 本章主要介绍一阶常微分方程初值问题的数值解法。
进一步: 令
y n1 y n
xn 1 xn
y n 1 y( x n 1 ) , y n y( x n )
f ( x , y( x ))dx h f ( x n , y n )
宽
9
高
实际上是矩形法
数值分析
第9章 常微分方程初值问题数值解法
(3)
用Taylor多项式近似并可估计误差
解决方法：有的可化为显格式，但有的不行 18
数值分析
第9章 常微分方程初值问题数值解法
与Euler法结合，形成迭代算法 ，对n 0，2， 1，
( yn0 )1 yn hf x n , yn ( k 1) h ( yn1 yn f x n , yn f x n1 , ynk )1 2
7
数值分析
第9章 常微分方程初值问题数值解法
建立数值解法的常用方法
建立微分方程数值解法,首先要将微分方程离散 化. 一般采用以下几种方法: (1) 用差商近似导数
dy yx yx x x dx x y
n 1 n n 1 n
n
,
n
进一步: 令
yn1 y( xn1 ) , yn y( xn )
由 x0 , y0 出发取解曲线 y y x 的切线（存在！），则斜率
常微分方程与数值解法
常微分方程与数值解法数学是自然界中最美丽的语言之一，常微分方程是数学中的一个重要分支。
常微分方程是研究随着时间推移而发生的连续变化的数学模型，是许多科学领域的数学基础，如物理学、天文学、生物学、化学、经济学等。
通过对微分方程的求解，我们可以预测未来的变化和趋势，制定相应的政策措施和科学研究方向。
一、常微分方程的基本概念常微分方程是包含未知函数及其导数的方程。
一般形式为dy/dx=f(x,y)，其中y为未知函数，x为自变量，f(x,y)是已知函数，称为方程的右端函数。
常微分方程可以分为初值问题和边值问题。
初值问题是指求解微分方程时需要给出一个特定的初值y(x)=y0，边值问题是指给出方程在一些点的值，而求出未知函数在整个区间上的值。
二、常微分方程的解法常微分方程有许多解法，例如分离变量法、齐次方程、全微分方程、一阶线性方程、变量分离法等。
其中，变量分离法是最基本和最重要的方法之一。
变量分离法的基本思想是将微分方程的未知函数y和自变量x分开，变成dy/g(y)=f(x)dx的形式，然后对两边进行积分。
三、数值解法的发展与应用数值解法是通过数值计算来求解微分方程的，它主要包括欧拉法、改进欧拉法、龙格-库塔法等。
欧拉法最简单、最基本，但精度较低，适用于解决一些简单的微分方程。
改进欧拉法和龙格-库塔法则精度更高，适用于解决较为复杂的微分方程。
数值解法在科学技术中的应用广泛，如气象学、环境保护、物理学、化学等。
以生态学为例，许多生态系统的动态变化可以用微分方程描述，如种群增长、捕食捕获、竞争关系等。
数值解法可以在一定程度上预测未来的生态状态，有助于制定相应的生态保护措施。
四、结论在现代科学技术中，微分方程和数值解法已经成为不可或缺的工具之一。
通过微分方程的求解，可以预测未来的变化和趋势，制定相应的政策措施和科学研究方向。
数值解法则更加精细和灵活，能够解决更为复杂的微分方程，广泛应用于各个领域。
因此，学习微分方程和数值解法，不仅是数学爱好者的追求，更是科学技术工作者不可或缺的技能。
数值分析(25) 常微分方程初值问题的
(9-1)
只要 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存 在唯一解。
3、 中点欧拉公式 /* midpoint formula */
中心差商近似导数
y( x 2 ) y( x0 ) y( x1 ) 2h
y( x2 ) y( x0 ) 2h f ( x1 , y( x1 ))
y n1 yn1 2hf (xn , y n ) n 1, 2,
h yn1 yn ( f ( xn , yn ) f ( xn 1 , yn hf ( xn , yn )) 2
可得
h yn 1 yn h (1 )( xn yn ) 1 2
yn 1 0.905 yn 0.095 xn 0.1
x0
x1
y1
x0 x1
过点 ( x0 , y0 ) ，以 f ( x0 , y0 ) 为切线斜率的
切线方程为
y y0 f ( x0 , y0 )( x x0 )
用 y1 y0 f ( x0 , y0 )( x1 x0 ) y0 hf ( x0 , y0 ) 近似代替 y ( x1 )
若数值公式的局部截断误差 En 1 O ( h p 1 ) ，则称此数值 公式是 p 阶的， p 是正整数。
欧拉法的局部截断误差,由Taylor展开: En 1 y( xn 1 ) yn 1
常微分方程的数值解与解析解
一、 常微分方程的解析解常微分方程的解析解也就是常微分方程的精确解，也称为常微分方程的符号解；一般可理解为求微分方程的通解或者特解的解析式或表达式；但只有少数的微分方程存在解析解。
在MA TLAB 中，由函数dsolve()求解常微分方程（组）的解析解，其具体格式如下： X=dsolve(‘方程1’，‘方程2’，…‘方程n ’，‘初始条件’，‘自变量’)函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。
例1：求解常微分方程1dy dx x y =+的MA TLAB 程序为：dsolve('Dy=1/(x+y)','x')，注意，系统缺省的自变量为t ，因此这里要把自变量写明。
结果为：-lambertw(-C1*exp(-x-1))-x-1其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X 。
例2：求解常微分方程2'''0yy y -=的MA TLAB 程序为：Y2=dsolve('y*D2y-Dy^2=0’,’x’) 结果为：Y2 =[ exp((x+C2)/C1)][ C2]我们看到有两个解，其中一个是常数。
例3：求常微分方程组253t tdx x y e dt dy x y e dt ⎧++=⎪⎪⎨⎪--=⎪⎩通解的MA TLAB 程序为：[X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t')例4：求常微分方程组020210cos ,224,0t t t dx dy x t x dt dt dx dy y e y dt dt =-=⎧+-==⎪⎪⎨⎪++==⎪⎩通解的MA TLAB 程序为：[X,Y]=dsolve('Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2','y(0)=0')二、 常微分方程的数值解在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。
《数值分析》第5讲：常微分方程数值解
1 2
y(3) ( xn )2 p2h3
o(h4 )
对照
y( xn1 )
yn
y'( xn )h
1 2
y''( xn )h2
第五章：常微分方程数值解
得
212 p2
1
1
可解得
yn1 yn h(1K1
K1 f ( xn , yn )
2
K
2
)
K
2
f ( xn p , yn
phK1 )
12
第五章：常微分方程数值解
§5.3 Lunge-Kutta方法
依据精度要求的待定系数法
x xn
x n p n1
1、二阶Lunge-Kutta方法（P114-P116)
令 xn p xn ph( xn, xn1), 0 p 1 yn1 yn h(1K1 2 K2 ) 加权平均斜率
K1 f ( xn , yn ) y(x) 在点 ( xn , 的yn斜) 率
例 P107
初值问题
y' y 2x / y
y(0)
1
0 x 1
补 充 ： 一 阶 线 性 方 程dy P( x) y Q( x)的 解
dx
常数变易法
y e p( x)dx
Q( x)e
p( x )dx
dx
C
非线性方程 dy P( x) y Q( x) yn Bernoulli型方程 dx
yn1 yn h(1K1 2 K 2 )
yn h[1 y'( xn ) 2 y'( xn p )]
yn
h1 y'( xn ) 2h
1 y'( xn ) 1
数值分析--第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 的步长。
今后如⽆特别说明，我们总假定步长为常量。
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解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乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
数值分析李庆扬第9章常微分方程初值问题数值解法讲义.
② 由 x0 , y0 f x0 , y0 切线 P0P1 ，
切线与 x x1 交点 P1 ： y1 的近似值 ；
③ 再由 x1 , y1 向前推进到 P2 ， 得到折线 P0P1 Pn ，近似 y yx 。
7
2021年5月4日
《数值分析》 黄龙主讲
h
yxn
yxn1
yn1 yn f
h
xn1 , yn1
yn1 yn h f xn1 , yn1
——后退的欧拉公式（隐式）
注意：① 显式计算方便，隐式稳定性较好；
② 上式隐含 yn1 ? ，采用迭代法求解。
12
2021年5月4日
《数值分析》 黄龙主讲来自欧拉公式的另一种理解：
将常微分方程 y f x, y 改写 dy f t , ytdt
“步进式”：顺着节点排列顺序，一步一步地向前推进。
步长：常用等步长 hn xn1 xn ，节点为 xn x0 nh 单步法：计算 yn1 时，只用到前一点的值 yn k 步法：计算 yn1 时，用到前面 k 点的值 yn , yn1 , , ynk1
5
2021年5月4日
《数值分析》 黄龙主讲
对微分方程从 xn 到 xn1 积分
y xn1 yxn
xn1 f t , yt dt
xn
由积分左矩形公式得
xn1 xn
f
t ,
yt dt
hf
xn ,
yxn
例如：
lim
h0
yxn1
h
yxn
yxn
yxn1
h
yxn
yxn
f xn , yxn
数值分析第九章常微分方程数值解法
数值分析第九章常微分方程数值解法常微分方程数值解法是数值分析中非常重要的一部分内容。
常微分方程是描述自然现象中动态变化规律的数学模型，解常微分方程可以揭示系统的变化趋势和规律。
然而，大多数常微分方程是无法通过解析方法求出解的，因此需要借助计算机进行数值计算。
数值解常微分方程方法主要包括：Euler方法、改进的Euler方法、四阶Runge-Kutta方法和龙格-库塔方法。
Euler方法是最简单的方法之一，它采用的是一阶Taylor展开式。
将待求的函数值与函数的一阶导数值代入Taylor展开式中，可以得到函数值在下一个时间步长上的近似值。
Euler方法的优点是简单易于实现，但其精度不够高，容易积累误差。
改进的Euler方法是对Euler方法的改进，它通过使用中间点上的导数值来减小误差。
改进的Euler方法的精度相比Euler方法要高一些，但仍然不够高。
四阶Runge-Kutta方法是目前使用较为广泛的数值解常微分方程的方法之一、它通过计算不同时间点上的斜率来估计函数值，在多个时间点上计算斜率的平均值来提高精度。
四阶Runge-Kutta方法的精度比Euler方法和改进的Euler方法要高，但计算量也相对较大。
龙格-库塔方法是数值解常微分方程中最常用的方法之一、它是四阶Runge-Kutta方法的延伸，通过计算不同时间点上的斜率来估计函数值，然后利用这些估计值计算更准确的斜率，在不同步长上进行迭代计算，直到满足所需精度。
龙格-库塔方法的精度比四阶Runge-Kutta方法要高，但计算量也相对较大。
除了以上几种方法外，还有一些其他数值解常微分方程的方法，如Adams法、Gear法等。
这些方法在不同场景下有着不同的适用性和优劣势。
总结起来，数值解常微分方程方法是研究常微分方程数值计算中的重要内容。
不同的方法有着不同的精度和计算量，可以根据具体问题的特点选择合适的方法进行数值计算。
然而，需要注意的是，数值解只是在给定的步长下对函数的近似值，可能会引入误差。
数值分析教案_常微分方程初值问题数值解法
第九章常微分方程初值问题数值解法图9-1n 作为()n x y 的近似值，得 ()n n y x hf ,)y x ,两边从n x 到1+n x 积分，得()dx x y x f x y x n nx x n n ⎰+=-+1))(,()1 矩形公式计算上式右侧积分，即()()x x x x x d x y x f dx x y x f n nn n⎰⎰++≈11,))(,()n ，得()n n n n y x hf y y ,1+=+，故欧拉法也称为矩形法。
为了达到较高精度的计算公式，对欧拉法进行改进，用梯形公式计算()()([1,2))(,(1++≈+n n n x f x y x f hdx x y x f n 的近似值，得9.2 龙格—库塔法前面讨论的欧拉法与改进的欧拉法都是一步法，即计算y 1+n 时，只用到前一步值。
龙格—库塔(Runge-Kutta)法(简称为R-K 方法)不是通过求导数的方法构造近似公式，而是通过计算不同点上的函数值，并对这些函数值作线性组合，构造近似公式，再把近似公式与解的泰勒展开式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。
我们先分析欧拉法与预估—校正法。
对于欧拉法⎩⎨⎧=+=+),(111n n n n y x hf k k y y 每步计算f 的值一次，其截断误差为O (2h )。
对于预估—校正法()()⎪⎪⎩⎪⎪⎨⎧++==++=+121211,,2121k y h x hf k y x hf k k k y y n n n n n n 每步计算f 的值两次，其截断误差为O (3h ).下面对预估—校正法进行改进，将该公式写成更一般的形式()()bh y ah x hf k y x hf k k R k R y y n n n n n n ++==++=+,,2122111 (2.1)其中b a R R ,,,21为待定常数。
选择这些常数的原则是在)(n n x y y =的前提下，使11)(++-n n y x y )的阶尽量高。
数学考研数值分析基础知识点
数学考研数值分析基础知识点数值分析是数学的一个分支，主要研究利用计算机进行数值计算的方法和算法。
在数学考研中，数值分析是一个重要的考点，本文将介绍数值分析的基础知识点，帮助考生能够更好地应对数值分析的考试。
一、插值与逼近1. 拉格朗日插值拉格朗日插值是一种常用的插值方法，用于根据已知数据点推测出未知点的值。
其基本思想是构造一个满足已知数据点的条件的拉格朗日多项式，并通过该多项式求解未知点的值。
2. 牛顿插值牛顿插值是另一种常用的插值方法，与拉格朗日插值相比具有更高的精度。
牛顿插值利用差商的概念，通过已知数据点的差商构造插值多项式，并利用该多项式求解未知点的值。
3. 最小二乘逼近最小二乘逼近是一种通过最小化残差平方和的方法，用于找到一个函数来近似已知数据点。
该方法常用于求解数据拟合问题，通过最小二乘逼近可以得到最优的拟合曲线。
二、数值积分1. 梯形公式梯形公式是一种常用的数值积分方法，通过将待积函数在积分区间上用一系列梯形逼近，从而求解积分的近似值。
梯形公式简单易懂，但精度比较低。
2. 辛普森公式辛普森公式是一种更高精度的数值积分方法，通过将待积函数在积分区间上用一系列二次曲线逼近，从而求解积分的近似值。
辛普森公式相比于梯形公式，在相同节点数的情况下有更高的精度。
三、常微分方程数值解法1. 欧拉法欧拉法是一种常用的常微分方程数值解法，通过将常微分方程转化为差分方程，从而近似求解方程的解。
欧拉法简单易懂，但对于某些情况下可能存在的数值不稳定性需要注意。
2. 修正的欧拉法和改进的欧拉法修正的欧拉法和改进的欧拉法是欧拉法的改进版，通过引入更高阶的项来提高精度和数值稳定性。
3. 4阶龙格-库塔法4阶龙格-库塔法是一种更高精度的常微分方程数值解法，通过迭代求解不同的插值点，并利用加权平均的方式来提高解的精度。
四、线性代数方程组的数值解法1. 直接法直接法是解线性代数方程组的一种常用方法，包括高斯消元法和LU分解法。
数值分析21(常微分方程数值解)
下面以2级方法为例子具体介绍龙格-库塔法。 yn1 yn hc1 K 1 hc2 K 2
K 1 f ( x n , yn ) K 2 f ( xn 2 h, yn h21 K 1 )
构造的基本思想是选择适当的系数使得方法的局部 截断误差阶数尽可能高。
Tn1 y( xn1 ) yn hc1 f ( xn , yn ) hc2 f ( xn 2h, yn h21 f ( xn , yn ))
y ( xn1 ) y ( xn ) h
f ( , y( )), 其中 [ xn , xn1 ]
h yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 2 是否可以推广改进的Euler方法？
(Runge-Kutta)龙格-库塔法
RK方法是一大类的方法, 其基本思想是采用如下形式:
yn1 yn hf ( xn1 , yn1 ( xn1 )) 隐式Euler法
梯形公式:
x n 1 xn
h f ( x , y( x ))dx [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
20:24
h yn1 yn [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
c1 c2 1, c22 1 / 2, c221 1 / 2
故方程组有无穷多组解, 每一组来自构成的RK公式的阶数 都是2, 都叫做二阶RK公式。
c1 c2 1 / 2, 2 21 1就是改进Euler法。
20:24
16/27
三阶Range-Kutta公式一般形式
《数值分析》 23
数值分析2.5__常微分方程的数值解
三、单步法的局部截断误差和精度 有关） 单步法的一般形式为 ： （ϕ与f 有关） 显式单步法 + hϕ ( x n , y n , y n +1 , h)
yn +1 = yn + hϕ ( xn , yn , h)
整体截断误差： 整体截断误差：从x0开始，考虑每一步产生的 开始， 误差，直到x 误差，直到 n，则有误差 e = y ( x ) − y
y′′ = f ( x, y, y′) a ≤ x ≤ b (3) y ( a ) = y0 , y (b) = yn
(1)，(2)式称为初值问题，(3)式称为边值问题。 ， 式称为初值问题 式称为初值问题， 式称为边值问题 式称为边值问题。 在实际应用中还经常需要求解常微分方程组： 在实际应用中还经常需要求解常微分方程组：
xn 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 y(xn) 0 0.1923 0.3448 0.4412 0.4878 0.5000 0.4918 0.4730 yn 0 0.2000 0.3840 0.5170 0.5824 0.5924 0.5705 0.5354 yn-y(xn) 0 0.0077 0.0392 0.0758 0.0946 0.0924 0.0787 0.0624
并将数值解和该问题的解析解比较。 并将数值解和该问题的解析解比较。 x 解析解： y ( x) = 2 1+ x 方法的具体格式： 解：Euler方法的具体格式： 方法的具体格式
y n +1
yn 2 = y n + h( − 2 y n ) xn
取h=0.2, xn=nh,(n=0,1,2…,15), f(x,y)=y/x – 2y2 计算中取f(0,0)=1. 计算结果如下： 计算结果如下： 计算中取
常微分方程初值问题的数值解法
常微分方程初值问题的数值解法在实际应用中，对于某些微分方程，我们并不能直接给出其解析解，需要通过数值方法来求得其近似解，以便更好地理解和掌握现象的本质。
常微分方程初值问题（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}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法，其需要计算的中间值较多，但是具有更高的精度。
实验八 常微分方程初值问题数值解法报告
实验八 常微分方程初值问题数值解法一、基本题科学计算中经常遇到微分方程（组）初值问题，需要利用Euler 法，改进Euler 法，Rung-Kutta 方法求其数值解，诸如以下问题：（1） ()⎪⎩⎪⎨⎧=-='004y xy y x y 20≤<x分别取h=0.1,0.2,0.4时数值解。
初值问题的精确解245x y e -=+。
（2） ()⎩⎨⎧=--='0122y y x y 01≤≤-x用r=3的Adams 显式和预 - 校式求解取步长h=0.1，用四阶标准R-K 方法求值。
（3）()()()100010321331221==-='⎪⎩⎪⎨⎧-='-='='y y y y y y y y y 10≤≤x用改进Euler 法或四阶标准R-K 方法求解取步长0.01,计算(0.05),(0.1y y y 数值解，参考结果 123(0.15)0.9880787,(0.15)0.1493359,(0.15)0.8613125y y y ≈-≈≈。
（4）利用四阶标准R- K 方法求二阶方程初值问题的数值解（I ）()()⎩⎨⎧='==+'-''10,00023y y y y y 02.0,10=≤≤h x(II)()()()⎩⎨⎧='==+'--''00,10011.02y y y y y y 1.0,10=≤≤h x(III)()()⎪⎩⎪⎨⎧='=+='00,101y y e y y x 1.0,20=≤≤h x(IV)()()⎩⎨⎧='==+''00,100sin y y y y 2.0,40=≤≤h x二、应用题1. 小型火箭初始质量为900千克，其中包括600千克燃料。
火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
许多实际问题的数学模型是微分方程或微分方程的定解问题。
如物体运动、电路振荡、化学反映及生物群体的变化等。
常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。
若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。
因此研究一阶微分方程的初值问题⎪⎩⎪⎨⎧=≤≤=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。
此方程称为差分方程。
需要说明的是，用不同的差商近似导数，将得到不同的计算公式。
(2) 用数值积分法将问题(9-1)中的微分方程在区间],[1+n n x x 上两边积分，可得)1,,1,0( ))(,()()(11-==-⎰++N n dx x y x f x y x y n nx x n n (9-3)用1+n y ,分别代替)(1+n x y ,)(n x y ，若对右端积分采用取左端点的矩形公式，即),())(,(1n n x x y x hf dx x y x f n n≈⎰+同样可得出显式公式(9-2)。
类似地，对右端积分采用其它数值积分方法，又可得到不同的计算公式。
(3) 用Taylor 多项式近似。
把1()n y x +在点处Taylor 展开，取一次多项式近似，则得2121()()()()2!()(,())() [,]2!n n n n n n n n h y x y x hy x y hy x hf x y x y x x ξξξ++'''=++''=++∈ 设1h ，略去余项，并以代替()n y x ，便得 1(,)n n n n y y hf x y +=+以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。
其中Taylor 展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。
上面我们给出了求解初值问题(9-1)的一种最简单的数值公式(9-2)。
虽然它的精度比较低，实践中很少采用，但它的导出过程能较清楚地说明构造数值解公式的基本思想，且几何意义明确，因此它在理论上仍占有一定的地位。
1简单的数值方法和基本概念1.1 Euler 法与向后Euler 法一、Euler 法Euler 方法就是用差分方程初值问题10(,) (0,1,,1)()n n n n y y hf x y n N y y a +=+=-⎧⎨=⎩(9-4)的解来近似微分方程初值问题(9-1)的解，即由公式(9-4)依次算出()n y x 的近似值(1,2,)n y n =。
从几何上看，微分方程(,)y f x y '=在xoy 平面上确定了一个向量场：点(,)x y 处的方向斜率为(,)f x y 。
问题(9-1)的解()y y x =代表一条过点00(,)x y 的曲线，称为积分曲线，且此曲线上每点的切向都与向量场在这点的方向一致。
从点000(,)P x y 出发，以00(,)f x y 为斜率作一直线段，与直线1x x =交于点111(,)P x y ，显然有1000(,)y y hf x y =+，再从出发，以11(,)f x y 为斜率作直线段推进到2x x =上一点222(,)P x y ，其余类推，这样得到解曲线的一条近似曲线，它就是折线012P PP 。
因此Euler 方法又称为Euler 折线法。
二、向后Euler 法在微分方程离散化时，用向后差商代替导数，即11()()()n n n y x y x y x h++-'≈，则得到如下差分方程11100(,) (0,1,,1)()n n n n y y hf x y n N y y x +++=+=-⎧⎨=⎩(9-5)用这组公式求问题(9-1)的数值解称为向后Euler 法。
向后Euler 法与Euler 法形式上相似，但实际计算时却复杂得多。
Euler 法计算1+n y 的公式中不含有1+n y ，这样的公式称为显式公式；向后Euler 法计算1+n y 的公式中含有1+n y ，称为隐式公式。
显式公式与隐式公式各有特点。
显式公式的优点是使用方便，计算简单，效率高。
其缺点是计算精度低，稳定性差；隐式公式正好与它相反，它具有计算精度高，稳定性好等优点，但求解过程很复杂，一般采用迭代法。
为了结合各自的优点，通常将显式公式与隐式公式配合使用，由显式公式提供迭代初值，再经隐式公式迭代校正。
上面隐式公式中，在求解1n y +时，1,n n y x +为已知，1n y +是方程111(,)n n n n y y hf x y +++=+的根。
一般说来，这是一个非线性方程，因此我们通过构造简单迭代法来求解。
迭代格式为(0)1(1)()111(,)(,) (0,1,2,)n n n n k k n n n n y y hf x y y y hf x y k +++++⎧=+⎪⎨=+=⎪⎩ 由于(,)f x y 满足Lipschitz 条件，所以(1)()()11111111(,)(,)k k k n n n n n n n n yy h f x y f x y hL yy +++++++++-=-≤-由此可知，只要1hL <，迭代法就收敛到解1n y +。
1.2 梯形公式利用数值积分方法将微分方程离散化时，若用梯形公式计算式(9-3)中右端积分，即111(,())[(,())(,())]2n nx n n n n x hf x y x dx f x y x f x y x +++≈+⎰ 并用1,n n y y +代替1(),()n n y x y x +，则得计算公式111[(,)(,)]2n n n n n n hy y f x y f x y +++=++(9-6)这就是求解初值问题(9-1)的梯形公式。
梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为(0)1(1)()111(,)(,)(,) (0,1,)2n n n n k k n n n n n n y y hf x y h y y f x y f x y k +++++⎧=+⎪⎨⎡⎤=++=⎪⎣⎦⎩ 由于函数(,)f x y 关于满足Lipschitz 条件，所以(1)()()11111111(,)(,)22k k k n n n n n n n n h hL y y f x y f x y y y +++++++++-=-≤- 其中为Lipschitz 常数。
因此，当012hL<<时，迭代法就收敛到解1n y +。
1.3 局部截断误差与方法的精度为了刻画近似解的准确程度，引入局部截断误差与方法精度的概念。
定义 假设在某一步的近似解是准确的，即()n n y y x =（这个假设称为局部化假设）。
在此前提下，用某公式推算所得1n y +，我们称111()n n n R y x y +++=-为该公式（即该方法）的局部截断误差。
定义9.2 如果某种方法的局部截断误差是1111()()p n n n R y x y O h++++=-=则称该方法是阶方法，或具有阶精度。
显然越大，方法的精度越高。
1）Euler 法的截断误差假设问题的解)(x y 充分光滑，且前步计算结果是精确的，即)(i i x y y =，()(,())i i i y x f x y x '=)(n i ≤于是Euler 法的截断误差是11111()()(,)()()()n n n n n n n n n n R y x y y x y hf x y y x y x hy x +++++=-=--'=--23()()2n h y x O h ''=+(9-7) 这里2()2n h y x ''称为局部截断误差主项。
显然21()n R O h +=。
2）向后Euler 法的截断误差。
计算公式是),(111++++=n n n n y x hf y y将(,)f x y 对用微分中值定理，有1111111(,)(,())(,)(())n n n n y n n n f x y f x y x f x y y x η+++++++=+-（在1n y +与1()n y x +之间）将11(,())n n f x y x ++在处Taylor 展开2111(,())()()()()n n n n n f x y x y x y x hy x O h +++''''==++于是231111()()()(,)(())()n n n n y n n n y y x hy x h y x hf x y y x O h η++++'''=+++-+将方程的解作Taylor 展开231()()()()()2n n n n h y x y x hy x y x O h +'''=+++因此2311111()()(,)(())()22n n n y n n n h hy x y y x f x y y x O h η+++++''-=---+故2311112311()()()1(,)21(,)()()2n n n n y n y n n h R y x y y x O h hf x h hf x y x O h ηη+++++⎡⎤''=-=-+⎢⎥-⎣⎦⎡⎤''⎡⎤=++-+⎢⎥⎣⎦⎣⎦23()()2n h y x O h ''=-+(9-8)3）梯形法的计算公式是111[(,)(,)]2n n n n n n hy y f x y f x y +++=++将(,)f x y 对用微分中值定理，有1111111(,)(,())(,)(())n n n n y n n n f x y f x y x f x y y x η+++++++=+-（在1n y +与1()n y x +之间）将11(,())n n f x y x ++在处Taylor 展开23111(,())()()()()()2n n n n n n h f x y x y x y x hy x y x O h +++'''''''==+++图9-3于是2341111()()()()(,)(())()242n n n n n y n n n h h hy y x hy x y x y x f x y y x O h η++++''''''=++++-+将方程的解作Taylor 展开2341()()()()()()26n n n n n h h y x y x hy x y x y x O h +''''''=++++因此3411111()()(,)(())()122n n n y n n n h hy x y y x f x y y x O h η+++++'''-=---+ 故341111341341()()()1(,)12[1(,)]()()12()()12n n n ny n y n n n h R y x y y x O h hf x h hf x y x O h h y x O h ηη+++++⎡⎤'''=-=-+⎢⎥-⎣⎦⎡⎤'''=++-+⎢⎥⎣⎦'''=-+(9-9)所以梯形法是二阶方法。