数值分析-常微分方程初值问题的解法
浅谈常微分方程初值问题数值解法
浅谈常微分方程初值问题数值解法在自然科学、工程技术、甚至社会科学的一些领域中，常常会遇见一阶常微分方程的求解问题：()上述问题，寻求解的具体表达式十分困难，仅对一些特殊形式的才有可能找到解的解析表达式，在大多情况下，初值问题的解不能用初等函数表示出来即使可写出解的解析表达式，但因为这些表达式过于复杂，要计算它在某些点上的函数值也异常困难。
在实际问题中，经常需要的恰是解在某些点上的函数值，因此研究初值问题的数值解法十分必要。
1 常微分方程初值问题的数值解法常微分方程的近似解法大体可分成三大类：一类是图解法和器械法；第二类是解的近似法；第三类是数值解法，即通过离散化的方法直接求出函数在某些点上的近似值，此数值解仅为精确解的近似解。
其基本原理为：一阶常微分方程的初值问题的解是上变量的连续函数，因此求上述问题的数值解，就是在区间上的若干离散点上用离散化的方法将初值问题化成离散变量的相应问题，从而相应问题的解可作为初值问题理论解的近似值。
由常微分方程的理论可知，只要在区域内连续，且关于满足林普希兹条件，则方程的解存在且唯一。
初值问题的数值解法通常采取“步进法”，而“步进法”又可分为“单步法”和“多步法”两类。
（1）单步法。
所谓“单步法”是指在计算时，只用到前一步的有关信息。
其一般形式为：，主要包括下面三种方法：Euler方法，改进的Euler公式-梯形公式和Runge-Kutta法。
（2）线性多步法。
单步法没有用到前几步计算得到的信息，因此为了提高精度，需重新计算多个点处的函数数值，如RK方法，故计算量较大。
线性多步法的基本思想是充分利用前面的已知信息来构造精度高且计算量小的算法来计算。
多步法常用方法是线性多步法，求解公式为：构造的常用方法是Taylor展开和数值积分方法。
常用的线性多步公式有：四阶Adams显式公式：四阶Adams隐式公式：四阶Milne显式公式：三阶Hamming公式：（隐式公式）预测校正系统和预测校正修正法：一般地，同阶的隐式法比显式法精确，而且数值稳定性好，但隐式公式中的求解较难，需要用到迭代法，这就增加了计算量。
第五章 常微分方程初值问题数值解法
则有
yn 1 yn hf ( xn , yn )
( 5.2 ) Euler格式
例5.1 用Euler格式解初值问题
2x y y y y (0) 1
取步长h=0.1.
(0 x 1)
Euler格式的具体形式为
y n 1 y n hf ( x n , y n ) 2 xn yn 0.1( yn ) yn 0.2 xn 1.1 yn yn
计算公式的精度 常以Taylor展开为工具来分析计算公式的精度. 为简化分析，假定yn是准确的，即在 yn y ( xn ) 的前提下估计误差 y ( xn 1 ) yn 1 Euler格式的局部截断误差 由 从而 局部截断误差
f ( xn , yn ) f ( xn , y ( xn )) y '( xn ) y ( xn 1 ) yn 1 y ( xn 1 ) ( yn hf ( xn , yn )) y ( xn 1 ) y ( xn ) hy '( xn )
y ( xn ), y ( xn 1 ), 的近似值 y1 , y2 , , yn , yn 1 ,
相邻两个节点的间距 h xi 1 xi 称为步长，步 长可以相等，也可以不等.本章总是假定h为定数， 称为定步长，这时节点可表示为
xn x0 nh , n 0,1, 2,
由f ( xn 1 , yn 1 ) f ( xn 1 , y ( xn 1 )) f y ( xn 1 , )( yn 1 y ( xn 1 )) f ( xn 1 , y ( xn 1 )) y '( xn 1 )（在xn点Taylor展开） h2 y '( xn ) hy ''( xn ) y '''( xn ) ... 2 3 2 h h 因此yn 1 y ( xn ) hy '( xn ) y ''( xn ) y '''( xn ) 2 4 h f y ( xn 1 , )( yn 1 y ( xn 1 )) 2 h2 h3 y ( xn 1 ) y ( xn ) hy '( xn ) y ''( xn ) y '''( xn ) 2 3!
数值分析 第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章 常微分方程初值问题数值解法
oa
b
a f ( x)dx (b a) f (b)
中矩形公式
b
ab
a f ( x)dx (b a) f ( 2 )
计算方法
梯形公式
bx
右矩形公式 中矩形公式 左矩形公式
§ 欧拉方法几何意义
y y y(x)
y0 y1 y2 0 x0 x1 x2
计算方法
x
§ 隐式欧拉方法
➢隐式欧拉法 /* implicit Euler method */
初 值 问 题 的 解 必 存 在 且唯 一 。
计算方法
§9.1 引言
三. 数值解法含义
所谓数值解法, 就是设法将常微分方程离散化, 建 立差分方程, 给出解在一些离散点上的近似值。
微分方程的数值解: 设方程问题的解y(x)的存在区 间是[a,b], 令a= x0< x1<…< xn =b, 其中hk=xk+1-xk, 如是等距节点h=(b-a)/n, h称为步长。
yi1 yi1 2h f ( xi , yi ) i 1, ... , n 1
计算方法
预估-校正法
三. 预估 — 校正法
/* predictor-corrector method */
方法 显式欧拉 隐式欧拉 梯形公式
中点公式
简单
稳定性最好
精度提高
精度低
精度低, 计算量大
计算量大
精度提高, 显式
在x0 x X上的数值解法。
四. 误差估计、收敛性
和稳定性
计算方法
§9.2 简单的数值方法与基本概念
一. 欧拉（Euler）格式
设 节 点 为xi a ih (i 0,1,2 , n) 方 法 一 ：Taylor展 开 法
常微分方程初值问题的数值解法
第七章 常微分方程初值问题的数值解法--------学习小结一、本章学习体会通过本章的学习，我了解了常微分方程初值问题的计算方法，对于解决那些很难求解出解析表达式的，甚至有解析表达式但是解不出具体的值的常微分方程非常有用。
在这一章里求解常微分方程的基本思想是将初值问题进行离散化，然后进行迭代求解。
在这里将初值问题离散化的方法有三种，分别是差商代替导数的方法、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 Ty 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 ε+=。
数值分析(25) 常微分方程初值问题的
忽略高阶项,取近似值可得到Euler公式
yn1 yn h f ( xn , yn ) (n 0,1, 2, ... )
数值分析
数值分析
3. 数值积分法区间 将方程y' f ( x, y)在区间 [ xn , xn1 ]上积分
xn1 y'dx xn1 f ( x, y)dx (n 0,1,L )
dy f ( x, y) x [a, b] dx y(a) y0
(9-1)
只要 f (x, y) 在[a, b] R1 上连续，且关于 y 满足 Lipschitz 条
件，即存在与 x, y 无关的常数 L 使 | f (x, y1) f (x, y2) | L | y1 y2 |
x1
记为
y1
过点 (x0 , y0 ) ，以 f (x0 , y0 ) 为切线斜率的 x0 x1
切线方程为 y y0 f (x0 , y0 )(x x0 )
用 y1 y0 f (x0 , y0 )(x1 x0 ) y0 hf (x0 , y0 ) 近似代替 y(x1)
数值分析
数值分析
h f (xn , yn )
y n 1
(n 0, 1, 2L )
数值分析
数值分析
例9-2 用改进的Euler方法解初值问题
y' y x 1
y(0)
1
取 h 0.1 ，计算到 x 0.5 。
解：利用
h yn1 yn 2 ( f (xn , yn ) f (xn1, yn hf (xn , yn ))
解：该问题的精确解为 y( x) y0e x
欧拉公式为 yn1 yn h yn (1 h) yn yn (1 h)n y0
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解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
数值分析 李庆扬 第9章 常微分方程初值问题数值解法
2017年1月4日
《数值分析》 黄龙主讲
讨论迭代的收敛性：
因函数
f x , y 对 y 满足利普希茨条件
f x , y1 f x , y2 L y1 y2 , y1 , y2 R
比较欧拉的后退公式和其 k 1 次迭代结果
yn1 yn h f xn1 , yn1
k 1 k yn y h f x , y 1 n n1 n1
两式相减得
k 1 k k yn y h f x , y f x , y hL y 1 n1 n1 n1 n1 n1 n1 yn1
由此可知：只要 hL 1 迭代法就收敛到解
yn1 时，只用到前一点的值 yn
k 步法：计算 yn1 时，用到前面 k 点的值 yn , yn1 , , ynk 1
5
2017年1月4日
《数值分析》 黄龙主讲
9.2
9.2.1
简单的数值方法
欧拉法与后退欧拉法
初值问题：
y f x , y
y x0 y0
解的形式：
由利普希茨条件，有
yn1 yn1
若选取 h 充分小，使得
20
k 1
hL k yn1 yn1 2
k hL 2 1 ，则 k 时有 yn1 yn1
2017年1月4日
《数值分析》 黄龙主讲
9.2.3
改进欧拉公式
① 先用向前欧拉公式，求得一个初步的近似值 预测：
称
f 关于 y 满足利普希茨条件， L 为 f 的利普希茨常数。
计算方法课件第八章常微分方程初值问题的数值解法
整体截断误差与局部截断误差的关系
定理：如果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方法
第一章常微分方程初值问题数值解法
（1.2.3）
其中rn,k(t)为插值余项。 代到（1.2.2）式中得
u ( tn +1 ) = u ( tn ) +
舍去余项 并用uj代替u(tj)即得
∫
tn+1 tn
Ln , k ( t ) dt + ∫ t rn , k ( t ) dt
tn+1
n
(1.2.4 (1.2.5)
Rn , k = ∫
⎡ ∑ ⎣α u
j =0
j =0
αk ≠ 0
（1.2.1）
j n+ j
⎤ − hβ j f n + j ⎦ = 0（数值解满足的差分方程）
因此称（1.2.1）为多步法 或 k-步法。 又因为(1.2.1)关于 u n + j , f n + j 是线性的，所以称为线性多步法。 为使多步法的计算能够进行，除给定的初值u0 外，还要 知道附加初值u1,u2,…,uk-1 ，这可用其它方法计算。 若 β k = 0 则称(1.2.1)是显式的； 若 β k ≠ 0 则方法（1.2.1）是隐式的。 例如，一般线性二步法可写成：
f ( t , u ( t ) ) = Ln , k +1 ( t ) + rn , k +1 ( t )
其中rn,k+1(t)为插值余项。 同理即
un +1 = un + h ∑ bk +1i f ( tn −i +1 , un −i +1 )
i =0
k +1
其中
bk +1i
=∫ ∏
−1
j =0 j ≠i
0
k +1
9、常微分方程初值问题数值解法
( k +1) yn +1
xn +1 ∫xn
− yn +1 |≤
hL 2
|
(k ) yn +1
− yn +1 |,
( 只要 hL < 1，则( 2.8)的ynk +1)收敛到(2.7)的yn +1. +1 2
三、单步法的局部截断误差与阶
一阶常微分方程初值问题(1.1)(1.2)的单步法的一般形式 yn +1 = yn + hϕ ( xn , yn , yn +1, h).
clear x=0,yn=1 %初始化 for n=1:10 yp=yn+0.1*(yn-2*x/yn); %预测 x=x+0.1; yc=yn+0.1*(yp-2*x/yp) ; yn=(yp+yc)/2 %校正 end
( 2.2)
作业： 作业：P381, 1, 2(1).
龙格—库塔 库塔(Runge-Kutta)法 §3 龙格 库塔 法
进一步 y ( xn +1 ) = y ( xn ) + ∫
xn +1 xn
f ( x, y ( x))dx,
(Байду номын сангаас.3)
∫
⇒ 其中
xn +1 xn
f ( x, y ( x))dx ≈ h ∑ ci f ( xn + λi h, y ( xn + λi h)).
yn +1 = yn + hϕ ( xn , yn , h),
i =1
r
(3.4) (3.5) 欧拉法r = 1, p = 1.改进
常微分方程初值问题的的数值解法
本章讨论常微分方程初值问题的数值解法
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 )
第7讲 常微分方程初值问题数值解法
yn1 yn hf ( xn1 , yn1 ) n 0,1, 2,, N 1
称上述公式为向后Euler 公式。 向后Euler 公式为隐式格式，需要利用迭代法求解
10
数值分析
第7讲 常微分方程初值问题数值解法
Euler方法的几何意义
③ 图形解
y
•龙格-库塔法
o
x
2
数值分析
第7讲 常微分方程初值问题数值解法
§1 引言 初值问题及其数值解的概念
dy f ( x , y ); a x b 一阶常微分方程初值问题： d x ( ) y( x ) y 0 0
常用的一些 解析解法：
分离变量法、变量代换、 常数变易法、Lapalace变换等
(3)
用Taylor多项式近似并可估计误差
h
2
y( x n 1 ) y( x n h) y( x n ) hy' ( x n )
y' ' ( ) h
2
2! y' ' ( x n )
y( x n ) hy' ( x n )
2!
进一步: 令
y n 1 y( x n 1 ) , y n y ( x n )
计算结果
（1）步长h=0.1的数值解比较表
x 0 0.1
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
精确解 1 1.0048
1.0187 1.0408 1.0703 1.1065 1.1488 1.1966 1.2493 1.3066 1.3679
向前欧拉 1 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}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法，其需要计算的中间值较多，但是具有更高的精度。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Taylor级数复习
y y( x)
x x0
则：
1 1 2 ( n ) n y ( x )( = y x )( + y x ) ( x x ) + y ( x ) ( x x ) + + y ( x ) ( x x ) + 0 0 0 0 0 0 0 2 ! n !
2 n h h ( n ) y ( x ) y ( x + h ) y ( x ) h y ( x ) y ( x ) y ( x ) i + 1 i i i i i 2 ! n !
2 n h ( h ) ( n ) y () x y ( x h ) y ( x ) h y ( x ) y ( x ) y ( x ) i 1 i i i i i 2 ! n !
x x i 0 , ... , n 1 ) 通常取节点间距 h 为步长，通常采用等 i i 1 i(求解,应用广泛,具
有应用价值。
Taylor级数复习
( x) 函数 y y 在
x 点作 x 0 Taylor级数展开：
1 2 yx ( )=yx ( 0) + y ( x ) ( x -x ) + y ( x ) ( x -x ) 0 0 0 0 2 ! 1 (n n + + y )( x ) ( x -x ) 0 0 + n !
2 n h h ( n + 1 ) y ( x ) y ( x h ) y ( x ) h y ( x ) y ( x ) y ( x ) i + 1 i i i i i 2 ! n !
2 n h ( h ) ( n + 1 ) y () x y ( x h ) y ( x ) h y ( x ) y ( x ) y ( x ) i 1 i i i i i 2 ! n !
y y hf ( x , y ) n 1 n n 1 n 1
是隐格式，要迭代求解
y n 1 y n 1
( k 1 ) ( 0 )
y hf ( x ,y n n 1 n 1 )
( k )
可以由向前差商公式求出
3、中心差商公式
y ( x ) y ( x ) h n 1 n 1 y ' ( x ) y ' ' ( ) n 1 n h 2 y y hf ( x , y ) n 1 n 1 n 1 n 1
计算方法
第六章 常微分方程初值 问题的数值解法
计算方法课程组
华中科技大学数学与统计学院
§6 常微分方程数值解法 §6.1 基本离散方法
§6.2 Runge-Kutta方法 §6.3 线性多步法
§6.4 收敛性与稳定性
§6.1 基本离散方法
考虑一阶常微分方程的初值问题 :
dy f (x , y) dx (a ) y 0 y x [a ,b ]
例如：
=x+y , x [01 ,] y (0 )1 y
x
其解析解为： y x 1 2 e
但是, 只有一些特殊类型的微分方程问题能够得到用解析 表达式表示的函数解，而大量的微分方程问题很难得到其解 析解。 因此，只能依赖于数值方法去获得微分方程的数值解。
y =e , x [0,1] 例如： y(0) 1
-x2
其解析解为：
y 1 e d t x [, 0 1 ]
t2 0
x
很难得到其解析解
例如：
=x+y , x [01 ,] y (0 )1 y
其解析解为
y x 1 2 ex
只有一些特殊类型的微分方程问题能够得到用解析表达式
表示的函数解，而大量的微分方程问题很难得到其解析解。 因此，只能依赖于数值方法去获得微分方程的数值解。 要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值 y y ( x ) ( i 1 , ... , n ) i i
6.1.2 Euler公式 利用等距分割，数值微分来代替导数项，建立差分格式。 b- a xi = i 称为局部截断误差。 m 1、向前差商公式
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 2 h y ( x ) y ( x ) hf ( x , y ( x )) y ' ' ( ) n 1 n n n n 2 所以，可以构造差分方程
显然，这个误差在 逐步计算过程中会 传播，积累。因此 还要估计这种积累
y y hf ( x , y ) n 1 n n n
定义: 在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考 虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差。 定义: 若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。 收敛性： 考察局部误差的传播和积累
2 h y ( x ) y ( x ) hf ( x , y ( x )) y ' ' ( ) n 1 n n n n 2 y y hf ( x , y ) n 1 n n n
2、向后差商公式
y ( x ) y ( x ) h n 1 n y ' ( x ) y ' ' ( ) n 1 n h 2 y ( x ) y ( x ) h n 1 n f ( x , y ( x )) y ' ' ( ) n 1 n 1 n h 2 2 h y ( x ) y ( x ) hf ( x , y ( x )) y ' ' ( ) n 1 n n 1 n 1 n 2