常微分方程数值解

总目录
本章目录
4.1
4.2
4.3
4.4
4.5
4.3 龙格-库塔方法
本课程主要研究其实际应用,，故直接给出各 类龙格-库塔公式。
1、二阶龙格-库塔
y n 1
yn
h 2
(k1
k2 )
k1 f (xn , yn )
k2 f (xn h, yn hk1 )
(4-11)
其中c1=1/2, c2=1/2, a=1, b=1
ch5
(4-18)
当h不大时，c可近似地看作常数。然后将步
长h对折，即取h/2为步长，从出发经过两步计算求
y断(x误n+差1)的为近c(似h/2值)5，，记于为是就yn(h有1/ 2，) 每一步计算的局部截
dy
y2
dx
,
0.0 x 0.4
y(0) 1
解：y0 1, h 0.1
用下面的迭代公式，对每个点迭代4次，k=1,2,3,4。
y
(0) n 1
yn
hyn2
y
( k 1) n 1
yn
h 2
[
y
2 n
(
y (k ) n 1
)
2
]
总目录
本章目录
4.1
4.2
4.3
4.2.4 梯形公式
4.1
4.2
4.3
4.4
4.5
4.1 微分方程在化工中的应用
另一个在化工中常见的微分方程是物料冷却过程 的数学模型，其模型可用下式表示：
dT dt
k (T
T0 )
（4-3）
在微分方程中我们称自变量函数只有一个的微分 方程为常微分方程，自变量函数个数为两个或两 个以上的微分方程为偏微分方程。给定微分方程 及其初始条件，称为初值问题；给定微分方程及 其边界条件，称为边值问题。
8
1131.43 2500 3868.57
9
1076.69 2500 3923.31
10
1019.76 2500 3980.24
11
960.546 2500 4039.45
12
898.968 2500 4101.03
n
xn
yn
zn
13 834.926 2500 4165.07
14 768.324 2500 4231.68
,
y
n
0.1 x 0.8
y(0) 1
k1
取步长h=0.2，计算公式为：kk
2 3
1
yn
0.2 6
k1
2k2
y
2 n
cos
xn
yn 0.1k1 2 cosxn
yn 0.1k2 2 cosxn
2k3
0.1 0.1
k
4
k4 yn 0.2k3 2 cosxn 0.2
n
xn
yn
总目录
本章目录
4.1
4.2
4.3
4.2.4 梯形公式
4.4
4.5
梯形公式：yn1
yn
h 2
(
f
(xn
,
yn
)
f (xn1, yn1 ))
(4-9)
梯形公式也是隐式格式。
用显式的欧拉公式和隐式的梯形公式给出的一 次预估-校正公式：
y
n1
yn
hf
(xn ,
yn )
yn1
yn
h[ f 2
(xn , yn )
1 3
hk1
k3
f
xn
2 3
h,
yn
1 3
hk1
hk
2
k4 f xn h, yn hk1 hk2 hk3
(4-17)
总目录
本章目录
4.1
4.2
4.3
4.4
4.5
4.3 龙格-库塔方法
实例
例4.3：用四阶龙格—库塔公式(4-16)求解下面初
值解问：题
dy dx
y2
cosx
4.4
4.5
该方程的精确解是
y
1 1
x
计算结果如表4-2所示。
n
xn
yn
1
0.1
1.1118
2
0.2
1.2520
3
0.3
1.4311
4
0.4
1.6763
表4-2 计算结果
y(xn )
1.1111 1.2500 1.4326 1.6667
实例
yn y(xn )
0.0007 0.0020 0.0095 0.0004
15 699.056 2500 4300.94
16 627.019 2500 4372.98
17
552.1
2500 4447.9
18 474.183 2500 4525.82
19 393.151 2500 4606.85
20 308.877 2500 4691.12
21 221.232 2500 4778.77
数值解yn。在计算中约定y(xn)表示常微分方程准
确解的值，yn表示y(xn)的近似值。
总目录
本章目录
4.1
4.2
4.3
4.4
4.5
4. 2 欧拉（Euler）公式
4.2.1 向前欧拉公式
4.2.2 向后欧拉公式
4.2.3 中心欧拉公式
4.2.4 梯形公式
总目录
本章目录
4.1
4.2
4.3
2k3
k4
k1 f xn , yn
(1)
k2
f
xn
1 2
h,
y
n
1 2
hk1
(2)
k3
f
xn
1 2 h, yn
1 2
hk
2
k4 f xn h, yn hk3
(4-16)
yn1
yn
h 8
k1
3k 2
3k3
k4
k1 f xn , yn
k 2
f
xn
1 ቤተ መጻሕፍቲ ባይዱ h, yn
22 130.081 2500 4869.92
23 35.2845 2500 4964.72
24 -63.3042 2500 5063.3
t
图4-3 三种初始值的温度变化曲线
表4-1
总目录
本章目录
4.1
4.2
4.3
4.4
4.5
4.2.1 向前欧拉公式
实例
从表4-1可以看到当自热引起物体温度升高的速度 小于散热引起温度下降的速度，物体的温度随时 间而逐渐减少：当自热引起物体温度升高的速度 与散热引起温度下降的速度平衡时，物体的温度 保持不变；当自热引起物体温度升高的速度大于 散热引起温度下降的速度，物体的温度随时间而 增长。在图4-3中L1，L2，L3分别表示初始值 3500，2500和1500的三条温度变化趋势曲线。
总目录
本章目录
4.1
4.2
4.3
4.4
4.5
4.1 微分方程在化工中的应用
在化工模拟中主要碰到的是常微分方程的初值问题：
y'(x) f (x, y)
y(a)
y0
,(a x b)
（4-5）
或
dy
dx
f (x, y) ,(a
x b)
y(a) y0
对于大多数常微分方程的初值问题, 只能计算它的 数值解。常微分方程初值问题的数值解就是求y(x) 在求解区间[a,b]上各个分点序列xn，n =1,2,…,m的
第4章 常微分方程数值解
4.1 微分方程在化工中的应用 4.2 欧拉（Euler）公式 4.3 龙格-库塔方法 4.4 常微分方程组的数值解法 4.5 程序示例及应用
总目录
4.1
4.2
4.3
4.4
4.5
4.1 微分方程在化工中的应用
微分方程在化工中应用的简单而又典型的例子是 套管式换热器的稳态温度分布。首先作以下假设：
4.2
4.3
4.4
4.5
4.2.1 向前欧拉公式
实例
解：wn1 wn h(0.04wn 100 ) 1.04wn 100, h 1 计算结w0果分见别表以4x-01=。1500，y0=2500，z0=3500代入。
w
5000 4500 4000 3500 3000 2500 2000 1500 1000 500
式：类似地，可得到计算y(xn+1)近似值yn+1的计算公
yn1 yn1 2hf ( xn , yn ) (4-8)
y算公公n的出式式值y计(14，才算-8再能)，称用求为为中得提中心y高n心+格1精的格式度值式算，。。出也因按y可2此公, 用y，式3,向…要(4后。-先8)欧y用，1可拉其需用公它要向式公知前计式道欧算计yn拉。-1,
4.4
4.5
4.2.1 向前欧拉公式
下式为计算近似值的向前欧拉公式：
yn1 yn hf ( xn , yn )
（4-6）
总目录
图4.1 欧拉折线法
图4-2 欧拉折线法几何示意图
本章目录
4.1
4.2
4.3
4.4
4.5
4.2.1 向前欧拉公式
实例
例4.1：假定某物体的温度w因自热而产生的热量可以
f
(xn1, yn1 )]
(4-10)
上式也称为改进的欧拉公式，它可合并成：
yn1
yn
h(f 2
(xn , yn )
f
( xn1 ,
yn
hf (xn ,
yn ))
总目录
本章目录
4.1
4.2
4.3
4.2.4 梯形公式
4.4
4.5
实例
例4.2：请用预估-校正
公式（改进的欧拉公 式）解右面初值问题：
总目录
本章目录
4.1
4.2
4.3
4.4
4.5
4.2.2 向后欧拉公式
向后欧拉公式：
yn1 yn hf (xn 1 , yn1 )
（4-6）
式(4-7)是yn+1的非线性方程，即隐式欧拉公式，用
迭代法求得yn+1。初始值
y
(0) n 1
由向前欧拉公式提供。
最简单的迭代公式为：
y y
(0) n1
y
(k 1)
n1
n
y
hf ( n hf
xn , (xn
yn )
1
,
y (k ) n1
, )
k
0,1,2,,
直到
y (k 1) n1
y (k ) n1
给定精度
h充分小时，以上迭代收敛。记( y) yn hf (xn1, y) ，则
' ( y) hf y (xn1, y) h充分小时，可保证 hf y (xn1, y) hL 1，
k1 k 2
f (xn , yn )
f
(xn
h 2
,
yn
h 2
k1 )
(4-13) (2)
k3 f (xn h, yn hk1 2hk2 )
yn
1
k1
f
yn
xn
h 4
k1
, yn
3k3
k 2
f
xn
1 3
h,
yn
1 3
hk1
(4-14)
k3
f
xn
2 3 h, yn
2 3
hk
yxn
yn yxn
1 2 3 4
总目录
0.2
1.24789
1.24792
0.4
1.63762
1.63778
0.6
2.29618
2.29696
0.8
3.53389
3.53802
表4-3 计算结果
本章目录
0.00003 0.00016 0.00078 0.00413
4.1
4.2
4.3
4.4
4.5
1、套管内侧为液体，其温度只随套管的长度改变 而改变，忽略温度的径向变化；套管环隙为蒸汽， 其温度在任何位置均为恒定值，可认为是饱和蒸 汽的温度。
2、忽略套管内侧流体的纵向热传导。
3、在整个套管长度方向上，总传热系数K不变。
总目录
本章目录
4.1
4.2
4.3
4.4
4.5
4.1 微分方程在化工中的应用
2
总目录
y
n1
k1
f
yn
h 9
2k1
xn , yn
3k 2
4k3
(3)
k 2
f
xn
1 2 h, yn
1 2
hk1
k3
f
xn
3 4
h,
y
n
3 4
hk
2
本章目录
(4-15)
4.1
4.2
4.3
4.4
4.5
4.3 龙格-库塔方法
3、四阶龙格—库塔公式
y n1
yn
h 6
k1
2k2
0 0
5 10 15 20 25
n
xn
yn
zn
1
1460
2500
3540
2
1418.4 2500 3581.6
3
1375.14 2500 3624.86
4
1330.14 2500 3669.86
5
1283.35 2500 3716.65
6
1234.68 2500 3765.32
7
1184.07 2500 3815.93
其中L为李普希兹条件。
总目录
本章目录
4.1
4.2
4.3
4.4
4.5
4.2.3 中心欧拉公式
y(x)的在x=x1处的中心差商式：y'(x1)
y(x2 ) y(x0 ) 2h
又y'(x1) f (x1, y(x1)) ，可得到y(x2)的近似值y2计算公式：
y2 y0 2hf ( x1 , y1 )
蒸汽入口
流体入口 ，u，t0
流体出口 ，u，tL
冷凝液出口
图4-1 套管式换热器温度分布示意图
流入的热量+传入的热量-流出的热量=0 （4-1）
r 2C put
K
2rdl (TW
t) r 2CPu(t
dl
dt ) dl
0
dt dl
2K u CPr
( TW
t)
（4-2） （4-3）
总目录
本章目录
使物体在每秒钟内以4%的速度增长，同时该物体由
于散热可使其温度在每秒种内下降100k，则物体温度
随时间变化的微分方程：dw 0.04w 100
（t以秒为单位)
dt
分别以初始温x(0)=1500k，y(0)=2500k，z(0)=3500k 用欧拉公式预测24秒后的物体温度趋势。
总目录
本章目录
4.1
4.3 龙格-库塔方法 步长的选择
下面以四阶龙格-库塔方法为例，说明如何自动选 择步长，使计算结果满足给定精度的要求。
格-库设塔从公节式点方x法n出经发过，一先步以计h算为得步y长(x，n+1利)的用近四似阶值龙，
记为
y (h) n 1
，由于公式的局部截断误差是y(h5)，故有
y(xn1 )
y (h) n1
或
y n 1
yn
hk2
k1 f (xn , yn )
(4-12)
k
2
f (xn
h 2 , yn
h 2 k1 )
其中c1=0,c2=1,a=1/2,b=1/2。
总目录
本章目录
4.1
4.2
4.3
4.4
4.5
4.3 龙格-库塔方法
2、三阶龙格-库塔公式
y
n
1
yn
h 6
(k1
4k2
k3 )
(1)