初值问题的数值解法

3 要求 Ri y( xi 1 ) yi 1 O(h ) ，则必须有：
1 1 2 1 , 2 p 2
这里有 3 个未知 数， 2 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到，p 1, 1 2 1 就是改进的欧拉法。
Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开
K 2 f ( xi ph, yi phK1 ) f ( xi , yi ) phf x ( xi , yi ) phK1 f y ( xi , yi ) O( h2 )
y( xi ) phy( xi ) O(h2 )
d f ( x, y) dx 首先希望能确定系数 1、2、p，使得到的算法格式有 2阶 dy 精度，即在 yi y( xi ) 的前提假设下，使得 f x ( x, y) f y ( x, y) dx Ri y( xi 1 ) yi 1 O(h3 ) f x ( x, y) f y ( x, y) f ( x, y) y( x )
y( x0 ) y0 yk 1 yk h f ( xk , yk 1 )
, k 0,1,...
隐式欧拉法的求解: 利用迭代的思路进行.
yi 1 yi hf ( xi , yi 1 )
变换为
y
( k 1) i 1
yi hf ( xi , y )
y i 1 K1 K2

1 1 y i h K 1 K 2 2 2 f ( xi , yi ) f ( xi h, yi hK 1 )
步长一定是一个h 吗？
§2 Runge-Kutta Method

例1 函数f ( t , y ) = t y 在区域D0 = {( t , y ) | 1 ≤ t ≤ 2, −3 ≤ y ≤ 4}
关于y满足Lipschitz条件，相应的Lipschitz常数可取为L = 2
3 存在性定理 定理1 设函数f ( t , y )在凸集D ⊂ R 2中有定义，若存在常数
(7.2.7)
称为显式Runge-Kutta(龙格-库塔 )方法，简称R-K方法，
其中正整数N 称为R-K方法的级，所有ci , ai , bij 都是待定 常数。
根据定义(7.2.7)，N 级R-K方法(7.9)的局部截断误差为
Rn+1 = y( t n+1 ) − y( t n ) − h∑ ci ki
dy 其斜率为 = f ( t0 , y0 ) dt ( t0 , y0 ) 由 点 斜 式 写 出 切线 方 程 dy y = y0 + ( t − t0 ) = y0 + ( t − t0 ) f ( t0 , y0 ) dt ( x0 , y )
0
等步长为h,则t1 - t0 = h, 可由切线算出 y1 ： 则 y1 = y0 + hf ( t0 , y0 ) 按此逐步计算y( tn ), 在tn +1处的值 ： yn+1 = yn + hf ( tn , yn ) y 注意： 这是“ 注意 ： 这是 “ 折 yN 线法” 而非“ 线法 ” 而非 “ 切 线法” 线法 ” 除第一个 点是曲线切线外， 点是曲线切线外 ， 其他点不是切线 y2 而是折线(如右 y1 y0 图所示)。 图所示 。
பைடு நூலகம்
则称数值解法(7.5)为显式方法。否则，称数值解法(7.3) 为隐式方法。

向与已给向量场在该点的方向一致。
而欧拉递推公式 (9.7) 所定义的点 ( xi , yi ) 可以看作 是通过准确初始点 ( x0 , y0 ) 的曲线多边形的顶点，而 每根连线的方向由它左边终点的向量场决定。
y
y y ( x)
o
x0
x1
x2
图9.2
xn
x
由于欧拉方法是用一条折线近似地代替解 y ( x) ， 所以欧拉方法也叫做折线法。
(9.23) , n 1
称为 r 阶泰勒( Taylor )方法。当 r 1 时就是欧拉方法的 公式，现取 r 2 ，则有
y0 y( x0 ) 2 h yi1 yi hf ( xi , yi ) ( f x( xi , yi ) f ( xi , yi ) f y( xi , yi )) i 0,1,2, , n 1 2
令
i1 yi1 yi*1
(9.17)
* i 1
* * y y 其中, i 1 是 yi 1 的近似值，若 i 1 不大，则 i 1 可以
作为 y( xi1 ) 的近似值， 若 i 1 很大，y 就不能用。
定义用一个数值方法，求解微分方程
y y
其中, 是一个复常数， 对给定步长 h 0 ， 在计算 yi 时 引入了误差 i 。 若这个误差在计算后面的 yi k (k 1, 2, ) 中所带来的误差按绝对值均不增加， 就说这个数值方 法对于这个步长 h 和复数 是绝对稳定的。
(9.26)
(9.26)式中的 yi 1 可以用欧拉方法得到，即得
(0) y i 1 yi hf ( xi , yi ) ( k 1) (k ) y y hf ( x , y i i 1 i 1 ) i 1

隐式方程通常用迭代法求解，而迭代过程的实质 是逐步显式化.
上页 下页
设用欧拉公式
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

积分法:
yk 1 yk h f ( xk , yk ) y ( x0 ) y0
积分项利用矩形公式计算
(1) y( xk 1 ) y( xk )
xk 1
xk
f (t , y(t ))dt
(★)

xk 1
xk
f (t , y(t ))dt h f ( xk , yk ) y( xk 1 ) y( xk ) h f ( xk , yk )
引言
初值问题的数值解法:求初值问题的解在一系列节点的值 y ( xn )的近似值 yn 的方法.本章数值解法的特点:都是采用“步进 式”,即求解过程顺着节点排列的次序一步步向前推进. 常微分方程初值问题: dy f ( x, y ), x [a, b] dx y ( x0 ) y0
替 f (x , y)关于 y 满足Lipschitz条件. 除了要保证(1)有唯一解外,还需保证微分方程本身是稳定的,即 (1)的解连续依赖于初始值和函数 f (x , y). 也就是说, 当初始值 y0 及函数 f (x , y)有微小变化时, 只能引起解的微小变化.
注: 如无特别说明,总假设(1)的解存在唯一且足够光滑. 在 f 连续有界, 则 f (x , y)对变量 y 可微的情形下, 若偏导数 y 可取L为
也称折线法 x
2. 梯形法
若采用梯形公式计算(★)中的积分项,则有 h y ( xk 1 ) y ( xk ) [ f ( xk , y ( xk )) f ( xk 1 , y ( xk 1 ))] 2 h yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk 1 )] 2 称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做 法: (0) (0) 先用Euler法求出初值 yk ,1 即 ,将其代入梯形公式 yk 1 yk h f ( xk , yk ) 的右端,使之转化为显式公式,即 h ( l 1) (l ) yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk (☆ ) 1 )] 2

4
当r=0时,称为单步法,这恰好就是欧拉方法.
3 当r =1时,称为二步法,其局部截断误差为 O h .
Adams线性多步法的特点:
r r , 1. 对任意的 j 1
r j 0
2. Adams线性多步法的截断误差的来源是数值积分的截断误差, 而数值积分的截断误差来源于插值误差,因此阿达姆斯方法的 的局部截断误差是 O hr 2 3. Adams线性多步法与Runge-Kutta法比较,R-K每提高一阶精 度至少需要多计算一个函数值;而线性多步法每提高一阶精 度只需多用一个已知的数据,所以从这个角度看计算量较小, 但线性多步法必须在单步法的基础上做. 所以在实际计算时 常常把两者相结合. 13
i
再对右端的积分应用数值积分(用函数值的线性组合来近似积分)
r yi 1 yi h j f xi j , yi j r j 0


12
由于 j 仅与r , j 有关, 而与f , i 无关，所以可以事先制成表。
r
当r =2时,称为三步法,其局部截断误差为 O h

i i
这个方法称二阶泰勒方法
6
2. 梯形方法和改进欧拉方法
xi 1 x i
h y x dx y xi 1 y xi [ f xi , yi f xi 1, yi 1 ] 2 h yi 1 yi [ f xi , yi f xi 1 , yi 1 ] 2
常微分方程初值问题解的存在唯一性定理:
以后我们总假定给出的方程都满足该定理的条件. 微分方程初值问题数值解法的特点： 先把方程离散化,即在区间[a,b ]中插入一些节点(通常采用等 距节点) ba

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展 开 法

第九章常微分方程初值问题的数值解法§1 相关知识§2 几种简单的数值方法§3 Runge–Kutta方法§4 单步法的进一步讨论§5 线性多步法西北工业大学理学院欧阳洁1上述定理称为一阶常微分方程初值问题解的适定性（存在性、惟一性与稳定性）定理。

对所讨论的一阶常微分方程初值问题，本章假设该问题是适定的，即解析解y(x)在区间[a,b]上是存在、惟一，且具有充分的光滑度。

因此f(x,y(x))也充分光滑。

西北工业大学理学院欧阳洁3西北工业大学理学院欧阳洁5常微分方程初值问题的数值解法分为：①单(一)步法：计算时，只用到和，即前一步的值。

1+n y n y n n x x ,1+显式单步法的一般形式为②多步法：计算时，除用到和以外，还用到和，即用到前k 步的值。

p n x −)1;1,2,1(>−=−k k p y p n L 1+n y n y n n x x ,1+对单步法与多步法，有显式与隐式方法之分。

显式、隐式多步法的一般形式类似。

隐式单步法的一般形式为),,(1h y x h y y n n n n ϕ+=+),,,,(111h y x y x h y y n n n n n n ++++=ϕ数值解法建立的过程：通过一定的离散化方法,将连续性问题的求解转化为有限个离散节点上解析解近似值的求解。

常用的离散化方法：Taylor 展开法；差商直接代替微商；数值积分法。

西北工业大学理学院欧阳洁16设一般的单步法为：显式公式隐式公式定义为某一数值方法在处的整体截断误差。

111()n n n e y x y +++=−1n x +五单步法的局部截断误差和阶整体截断误差不仅与这一步的计算有关，还依赖于前面各步的计算。

1+n x 011,,,,x x x x n n L −下面着重分析计算中某一步的误差（局部截断误差），后面将对显式单步法，给出整体截断误差与局部截断误差之间的关系。

初值问题与解方法初值问题是数学中的一个重要概念，它涉及到微分方程的解的初始条件。

解决初值问题的方法有多种，本文将介绍几种常用的解法，并讨论它们的适用性和优缺点。

一、欧拉法（Euler's method）欧拉法是一种较为简单的数值解法，通过逐步逼近微分方程的解。

它的基本思想是将时间和空间分割成小的步长，并用线性逼近的方式求解微分方程。

欧拉法的计算公式为：y_{n+1} = y_n + h * f(t_n, y_n)其中，y_{n+1} 是下一个时间步长上的解，y_n 是当前时间步长上的解，h 是步长（时间或空间），f(t_n, y_n) 是微分方程的右端函数。

欧拉法的优点是简单易懂、计算量小。

然而，它的精度较低，对于具有较大步长或非线性的微分方程，可能会产生较大的误差。

二、改进的欧拉法（Improved Euler's method）改进的欧拉法是对欧拉法的一种改进方法，通过增加一个中间点的计算来提高精度。

改进的欧拉法的计算公式为：y_{n+1} = y_n + (h/2) * (f(t_n, y_n) + f(t_{n+1}, y_n + h * f(t_n, y_n)))改进的欧拉法通过使用两个不同的斜率来进行计算，提高了解的逼近精度。

相比于欧拉法，改进的欧拉法的精度更高，误差较小。

三、龙格-库塔方法（Runge-Kutta methods）龙格-库塔方法是一类常用的数值解法，包括一阶、二阶、四阶等不同精度的方法。

其中，最常用的是四阶龙格-库塔方法。

四阶龙格-库塔方法的计算公式为：k_1 = h * f(t_n, y_n)k_2 = h * f(t_n + h/2, y_n + k_1/2)k_3 = h * f(t_n + h/2, y_n + k_2/2)k_4 = h * f(t_n + h, y_n + k_3)y_{n+1} = y_n + (k_1 + 2k_2 + 2k_3 + k_4)/6四阶龙格-库塔方法通过使用多个斜率进行逼近，进而提高了解的精度。

常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解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乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。

整体截断误差与局部截断误差的关系
定理：如果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方法

初值问题数值解法
计算方法
第27讲 • 2018.12.19
– 1 / 23
9.1 初 值 问 题 常 用 数 值 方 法 的 建 立 与 使 用
⊲ 初值问题数值解法
基本数值方法的建立与隐 式法的求解 龙格- 库塔法 待定系数法、预测- 校正 公式
第 9章
常微分方程数值解法
9.1.1 基本数值方法的建立与隐式法的求解 仍然考虑一阶常微分方程初值问题 y ′ (x) = f (x, y (x)),a ≤ x ≤ b, y ( a ) = y0 . 无论通过数值微分法、数值积分法还是泰勒级数法, 所得到的数值算 法从求解方式上可分为显式法和隐式法, 从当前数值解与求解历史的 关系上又分为多步法和单步法. 对于多步法的实现, 涉及到开始步的 计算, 对于隐式法涉及到代数方程(尤其是非线性代数方程)的求解方 法. 下面我们就来考虑这些问题.
由y ′ (x) = f (x, y (x))知
′ ′ ′ ′ ′ y ′′ (x) = fx + fy y = fx + fy f, ′′ ′′ ′′ 2 ′ ′ ′2 y ′′′ (x) = fxx + 2fxy f + fyy f. f + fx fy + fy
第27讲 • 2018.12.19
当|yi+1
(k+1)
− yi+1 | ≤ ε时，取yi+1 ≈ yi+1 .
(ห้องสมุดไป่ตู้ )
第27讲 • 2018.12.19
– 5 / 23
9.1 初 值 问 题 常 用 数 值 方 法 的 建 立 与 使 用
初值问题数值解法 基本数值方法的建立 与隐式法的求解 龙格- 库塔法 待定系数法、预测- 校正 公式

第九章常微分方程初值问题的数值解法第一部分内容提要一、数值解的一般概念常微分方程初值问题00'()(,)()y x f x y y x y =⎧⎨=⎩的数值解是指通过一定的近似方法得出准确解()y y x =在一列离散点012,,,,,n x x x x 上的近似值012,,,,,n y y y y 。

数值解的特征是步进式，即()y x 在1n x +点的近似值1n y +是由1,,n n x x -等若干点处的近似值1,,n n y y -的信息给出的递推公式。

若1n y +依赖于前面k 步的值11,,,n n n k y y y --+，则称为k 步法；1k =称为单步法。

利用()y x 在11,,,n n n k x x x --+的精确解11(),(),,()n n n k y x y x y x --+借助某种算法计算出1n y +，则称11()n n y x y ++-为该方法的局部截断误差。

如果一个算法的局部截断误差是1()p O h +，则称该方法是p 阶的；而利用数值解11,,,n n n k y y y --+得到的1n y +与微分方程的精确解之差11()n n y x y ++-称为整体截断误差，即是该数值方法的误差。

对于固定的0x x >，取0x x h n-=，用某种算法得到n y ，如有lim ()n h y x y →-=0，则称该方法是收敛的。

注意，因x 是固定的，随着0h →，数值解的步数n →∞。

二、在实际计算时因为舍入误差不可避免，实际得到数值解是n y ，稳定性即研究n n y y -是否随着计算步骤n 的增加而增加。

通常所提的稳定性是通过模型方程(0)y y λλ'=<来讨论的。

若当某一步n y 有舍入误差时，在以后的计算中误差不会逐步扩大，则称这种稳定性为绝对稳定性。

三、简单单步法及其收敛性、稳定性Euler 法1(,)n n n n y y hf x y +=+的局部截断误差为2()O h ，整体截断误差为()O h ，即一阶收敛。

第7章 常微分方程初值问题的数值解法7.1 引 言科学研究和工程技术中的许多问题在数学上往往归结为微分方程的求解问题。

为了确定微分方程的解，一般要加上定解条件，根据不同的情况，这些定解条件主要有初始条件(initial condition )和边界条件(boundary condition). 只含初始条件作为定解条件的微分方程求解问题称为初值问题(initial-value problem); 例如天文学中研究星体运动,空间技术中研究物体飞行等，都需要求解常微分方程初值问题(initial-value problem for ordinary differential equations). 只含边界条件作为定解条件的微分方程求解问题称为边值问题(boundary-value problem).除特殊情形外，微分方程一般求不出解析解，即使有的能求出解析解，其函数表示式也比较复杂，计算量比较大，而且实际问题往往只要求在某一时刻解的函数值. 为了解决这个问题，有两种方法可以逼近原方程的解。

第一种方法是：将原微分方程化简为可以准确求解的微分方程，然后使用化简后的方程的解近似原方程的解；第二种方法是：将求原微分方程的解析解转化为求原方程的数值解，这是实际中最常用的方法。

本章将介绍求解常微分方程初值问题的常用的数值方法。

第8章将介绍常微分方程的边值问题的常用的数值方法。

为简明起见，本章主要介绍形如()(,),,()y t f t y a t b y a α'=≤≤⎧⎨=⎩(7.1.1)的初值问题的数值解法. 在介绍这些方法之前，需要了解常微分方程的一些相关定义和结果.定义7.1.1 函数(,)f t y 称为在集合2D ⊂R 上关于变量y 满足Lipschitz 条件(Lipschitz condition)，简称Lip 条件，如果存在常数0L >，使得1212(,)(,)f t y f t y L y y -≤- (7.1.2)对所有12(,),(,)t y t y D ∈都成立. 常数L 称为Lipschitz 常数(Lipschitz constant)，简称Lip 常数.定义7.1.2 如果对所有1122(,),(,)t y t y D ∈，都有()1122(1),(1))t y t y D λλλλ-+-+∈， (7.1.3)其中01λ≤≤，那么称集合2D ⊂R 为凸集(convex set).()y t .定理7.1.2 设函数(,)f t y2上. 如果存在常数0L >，使得(7.1.4)对一切(,)t y D ∈成立，那么(,)f t y 在D 上关于变量y 满足Lipschitz 条件，且L 为求微分方程数值解的主要问题：(1) 如何将原微分方程离散化, 并建立求其数值解的递推公式;(2) 如何求递推公式的局部截断误差, 数值解n y 与精确解()n y t 的误差估计； (3) 研究递推公式的稳定性与收敛性.7.2 Euler 方法及改进的Euler 方法7.2.1 Euler 格式与梯形格式考虑一阶常微分方程的初值问题()(,),,().y t f t y a t b y a α'=≤≤⎧⎨=⎩(7.2.1) 设011N N a t t t t b -=<<<= ，其中0,0,1,,n t t n h n N =+= 为等距节点，步长()h b a N =-.在1[,](0,1,,1)n n t t n N +=- 上对(,())y f t y t '=两边积分，得11()()(,())d n nt n n t y t y t f t y t t ++=+⎰(7.2.2)7.2.1.1 Euler 格式用左矩形求积公式计算式(7.2.2)右端积分项，得1(,())d (,())n nt n n t f t y t t h f t y t +≈⎰,代入式(7.2.2)右端，得1()()(,())n n n n y t y t h f t y t +≈+. (7.2.3)用()n y t 的近似值n y 代入式(7.2.3)右端，记所得结果为1+n y ，得10(,),0,1,,1,,n n n n y y h f t y n N y α+=+=-⎧⎨=⎩ (7.2.4) 并称式(7.2.4)为求解初值问题(7.2.1)的Euler 方法(Euler ’s method)或Euler 格式(Euler scheme)，1(,)n n n n y y h f t y +=+称为差分方程(difference equation).注 Euler 方法是最早的解决一阶常微分方程初值问题的一种数值方法，虽然它的精度不高，很少被采用，但是它反映了微分方程的数值解法的基本思想和特征.若式(7.2.2)右边的积分由数值积分的右矩形公式近似，并用近似值n y 替代()n y t ，近似值1+n y 替代1()n y t +，则可得到1110(,),0,1,,1,,n n n n y y h f t y n N y α+++=+=-⎧⎨=⎩ (7.2.5) 并称式(7.2.5)为后退的Euler 方法(backward Euler ’s method) 或后退的Euler 格式(backward Euler ’s scheme). 111(,)n n n n y y h f t y +++=+是差分方程。

本章讨论常微分方程初值问题的数值解法
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 )

然后用 y n 代替 y ( x n ) ，即得
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}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法，其需要计算的中间值较多，但是具有更高的精度。