第六章 6.1 常微分方程及其求解概述 6.2 初值问题解法 6.3 边值问题解法 第六章 常微分方程及 方程组的解法 6.1 常微分方程及其求解概述 6.2 初值问题解法 6.3 边值问题解法 浙江大学研究生学位课程 《实用数值计算方法》 常微分方程及方程组的解法
6.1 常微分方程及其求解概述 6.2 初值问题解法 6.3 边值问题解法 1.Euler方法 2.线性多步法 6.1 常微分方程及其求解概述 6.2 初值问题解法 1.Euler方法 2.线性多步法 3.Runge--Kutta法 4.方程组及刚性问题的Gear方法 6.3 边值问题解法 1.Shooting(试射法) 2.差分法 浙江大学研究生学位课程 《实用数值计算方法》
6.1 常微分方程及求解概述 （Ordinary Differential Equations, ODE) 6.1.1 基本概念 描述自由落体的ODE: 浙江大学研究生学位课程 《实用数值计算方法》
(6-2),(6-3)是(6-1)的初始条件。亦称定解 条件。(6-1)(6-2)叫做初值问题。 6.1.1 只有一个自变量的微分方程为 ODE, 否则称为偏微分方程 PDE。 方程中未知函数导数的最高阶数称为 方程的阶。 （6-4）是二阶的 方程中关于未知函数及其各阶导数 均是一次的，则称为线性微分方程。 和（6-1）都是线性二阶ODE。 (6-2),(6-3)是(6-1)的初始条件。亦称定解 条件。(6-1)(6-2)叫做初值问题。 浙江大学研究生学位课程 《实用数值计算方法》
（6-1），（6-3）叫做边值问题。 在没有给定解条件时。方程一般 有一族解曲线 y(x,c) 。如： 6.1.1 （6-1），（6-3）叫做边值问题。 在没有给定解条件时。方程一般 有一族解曲线 y(x,c) 。如： 对任意的n阶ODE，如果能写成： 则称该方程为显式的。方程（6-4）是 显式的。而下面方程是隐式的。 浙江大学研究生学位课程 《实用数值计算方法》
6.1.1 对于高阶显式方程。通过定义n-1个 新变量，可以写成n维一阶方程组。 即令： 浙江大学研究生学位课程 《实用数值计算方法》
当 f(x,y)与y无关时，f(x,y)=g(x) 6.1.1 在讨论初值问题时，我们从一阶方 程开始： 然后毫不费力地套用来解方程组。 当 f(x,y)与y无关时，f(x,y)=g(x) 浙江大学研究生学位课程 《实用数值计算方法》
6.1.2 数值解及其重要性 浙江大学研究生学位课程 《实用数值计算方法》
6.1.3 ODE数值解的基本思想和方法特点 基本思想有两点 1. 离散化 用Taylor级数，数值积分和差商 的代数方程(称差分方程)。 2. 递推化 在具有唯一解的条件下，通过 步进法逐步计算出解在一系列离散 点上的值。从而得到原ODE的数值 近似解。 浙江大学研究生学位课程 《实用数值计算方法》
6.2 初值问题解法 我们讨论一阶ODE，而高阶可 能化为一阶ODEs。一阶初值问题 可以一般地写成： 6.2.1 欧拉（Euler）方法 6.2 初值问题解法 我们讨论一阶ODE，而高阶可 能化为一阶ODEs。一阶初值问题 可以一般地写成： 6.2.1 欧拉（Euler）方法 Euler方法是求解(6-8)最简单方法， 但精度差，故不实用。然而对理论分 析很有用。 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 6.2.1.1 方法原理及推导 设初值问题(6-8)满足： 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 图 6.1 常微分方程初值问题的数值解 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 欧拉方法的几何意义： h步长 图 6.2 Euler方法的几何意义 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 单步法 自动起步 显式 多 步 法 隐式 半隐式 图 6.3 ODE求解方法的类型 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 表 6.1 自由落体运动方程的Euler公式求解 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 图 6.4 运动轨迹 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1.2 Euler方法的误差估计 一般其它方法的误差估计也类似。 这里误差是指截断误差(算法理论误 差) 而不是舍入误差。后者由计算机字 长等决定，属于稳定性问题。 i) 几何分析 图 6.5 Euler公式的误差 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 局部截断 误差， 浙江大学研究生学位课程 《实用数值计算方法》
整体截断误差：设ym是Euler公式（6-9） 精确解，而 y(x) 是初值问题（6-8）的解。 则 整体截断误差定义为 6.2.1 整体截断误差：设ym是Euler公式（6-9） 精确解，而 y(x) 是初值问题（6-8）的解。 则 整体截断误差定义为 它是局部截断误差的积累。 定理：若f(x,y)关于y满足Lipschitz条件。 则有估计式： 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 注意： 稳定性： 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 表 6.2 予估校正求解结果对比 浙江大学研究生学位课程 《实用数值计算方法》
6.2.1 表 6.3 Euler法与外推结果的比较 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 线性多步法 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 () () 浙江大学研究生学位课程 《实用数值计算方法》
表 6.4 外插系数bki值 Adams 外插法 (k=2) 3阶3步 显式 图 6.6 3阶3步外插法 6.2.2 浙江大学研究生学位课程 图 6.6 3阶3步外插法 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 Adams 外插法 (k=2) 4阶3步 图 6.7 4步3阶Adams内插公式 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 图 6.8 一般化插值形式 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 图 6.9 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 线性多步法的绝对稳定性： 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 定义： 绝对稳定。 绝对稳定区域。 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 Milne 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 表 6.6 计算结果 浙江大学研究生学位课程 《实用数值计算方法》
（Predictor--Corrector Method） 6.2.2 6.2.2.4 预估--校正方法 （Predictor--Corrector Method） 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 浙江大学研究生学位课程 《实用数值计算方法》
注意：一步校正的计算量 预估计算量。 所以要适当选取h才能发挥PC的优点。 设绝对稳定区域： 达到精度的校正次数为N，则h的选取， 6.2.2 注意：一步校正的计算量 预估计算量。 所以要适当选取h才能发挥PC的优点。 设绝对稳定区域： 达到精度的校正次数为N，则h的选取， 应满足： 否则可用N步显式算法稳定达到目的。 h 浙江大学研究生学位课程 《实用数值计算方法》
6.2.2 Adams四步四阶预估校正算法 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 Runge--Kutta 方法 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 六个未知数，二个自由，故可取 故： 浙江大学研究生学位课程 《实用数值计算方法》
----最常用的古典Runge-Kutta方法。 增加计算函数值的次数（级）与提高精 度（阶）的关系见下表： 6.2.3 N=4: 四级四阶R--K方法 ----最常用的古典Runge-Kutta方法。 增加计算函数值的次数（级）与提高精 度（阶）的关系见下表： 表 6.7 Runge-Kutta方法中级与阶的关系 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 表 6.8 各种解法在例题中的结果比较 浙江大学研究生学位课程 《实用数值计算方法》
（2）. 单步法，自动起步 （3）. 易改为变步长 （4）. 绝对稳定区域较同阶线性多步法大 （5）. 计算工作量较大，有时大于隐式方法 6.2.3 （2）. 单步法，自动起步 （3）. 易改为变步长 （4）. 绝对稳定区域较同阶线性多步法大 （5）. 计算工作量较大，有时大于隐式方法 （6）. 估计误差不易 绝对稳定性讨论 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 表 6.9 各级R-K方法的绝对稳定区域 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 表 6.10 不同步长对精度的影响 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 浙江大学研究生学位课程 《实用数值计算方法》
图 6.10 变步长Runge--Kutta方法框图 6.2.3 P阶 图 6.10 变步长Runge--Kutta方法框图 浙江大学研究生学位课程 《实用数值计算方法》
6.2.3 浙江大学研究生学位课程 《实用数值计算方法》
6.2.4 出发值的计算 浙江大学研究生学位课程 《实用数值计算方法》
6.2.4 浙江大学研究生学位课程 《实用数值计算方法》
质量控制 Runge--Kutta 步进程序 6.2.4 质量控制 Runge--Kutta 步进程序 SUBROUTINE RKQC(Y, DYDX, N, X, HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS) PARAMETER (NMAX=10, PGROW=-0.20, PSHRINK=-0.25, FCOR=1./15. ONE=1.0, SAFETY=0.9, ERRCON=6.E-4 EXTERNAL DERIVS DIMENSION Y(N), DYSX(N), YSCAL(N), YTEMP(NMAX), YSAV(NMAX), DYSAV(NMAX) XSAV=X 保留初值 DO 11 I=1，N YSAV(I)=Y(I) DYSAV(I)=DYDX(I) H=HTRY HH=0.5H CALL RK4 (YSAV, DYSAV, N, XSAV, HH, YTEMP, DERIVS) X=XSAV+HH CALL DERIVS (X, YTEMP, DYDX) CALL RK4 (YTMP, DYDX, N, X, HH, Y, DERIVS) X=XSAV+H IF (X. EQ. XSAV) PAUS E ‘步长无意义’ ERRMAX=0. DO 12 I=1，N YTEMP(I)=Y(I)-YTEMP(I) ERRMAX = MAX(ERRMAX, ABS(YTEMP(I) / YSCAL(I))) ERRMAX = ERRMAX / EPS 置初始值 11 1 两个半步长计算 一个单步计算 12 误差计算 浙江大学研究生学位课程 《实用数值计算方法》
COMMON /PATH/ KMAX，KOUNT，DXSAV，X(200)，Y(10，200) EXTERNAL DERIVS RKQC 6.2.4 PARAMETER (NVAR=4) DIMENSION VS(NVAR) COMMON /PATH/ KMAX，KOUNT，DXSAV，X(200)，Y(10，200) EXTERNAL DERIVS RKQC X1=1.0 X2=10.0 VS(1)=BESSO(X1) VS(2)=BESSI(X1) VS(3)=BESSJ(2,X1) VS(4)=BESSJ(3,X1) EPS=1.0E-4 HI=1.0 HMIN=0.0 KMAX=100 DXSAV=(X1-X2)/20.0 CALL ODEINT (VS，NVAR，X1，X2，EPS，HI，HMIN，NOK NBAD， NOK，NBAD，DERIVS，RKQC) WRITE ( ， ) END SUBROUTINE DERIVS (X，Y，DYDX) DIMENSION Y(1)，DYDX(1) DYDX(1) = -Y(2) DYDX(2) = Y(1)-(1.0/X) Y(2) DYDX(3) = Y(2)-(2.0/X) Y(3) DYDX(4) = Y(3)-(3.0/X) Y(4) RETURN 浙江大学研究生学位课程 《实用数值计算方法》
6.2.4 IF (ERRMAX.GT.ONE) THEN 不满足要求 H = SAFETY H (ERRMAX PSHRINK) 估计下次步长（缩小） GOTO 1 ELSE 满足要求 HDID = H IF (ERRMAX. GT. ERRCON) THEN HNEXT = SAFETY H (ERRMAX PGROW) 估计可放大的步长 ELSE HNEXT =4. H 最多可放大四倍 ENDIF DO 13 I=1，N 完成五阶截断误差 Y(I)=Y(I) + YTEMP(I) FCOR CONTINUE RETURN END 四阶Runge--Kutta 步进程序（续） 13 浙江大学研究生学位课程 《实用数值计算方法》
6.2.4 计算结果 NOK=30 NBAD=0 KOUNT=15 浙江大学研究生学位课程 《实用数值计算方法》
6.2.4 四阶 Runge--Kutta 方法： SUBROUTINE RK4(Y，DYDX，N，X，H，YOUT，DERIVS) PARAMETER (NMAX=10) DIMENSION Y(N)，DYDX(N)，YOUT(N)，YT(NMAX) DYT(NMAX)，DYM(NMAX) HH=H 0.5 H6=H/6.0 XH=X+HH DO 11 I=1，N YT(I)=Y(I)+HH DYDX(I) 第一步 CALL DERIVS(XH，YT，DYT) DO 12 I=1，N 第二步 YT(I)=Y(I)+HH DYT(I) CALL DERIVS(XH，YT，DYM) DO 13 I=1，N 第三步 YT(I)=Y(I)+H DYM(I) DYM(I)=DYT(I) + DYM(I) CALL DERIVS(X+H，YT，DYT) DO 14 I=1，N YOUT(I)=DYDX(I)+DYT(I)+2. DYM(I) YOUT(I)=Y(I) + H6 YOUT(I) CONTINUE RETURN END 11 12 13 14 浙江大学研究生学位课程 《实用数值计算方法》
用四阶 RUNGE--KUTTA公式，等步长计算NSTEP步 6.2.4 用四阶 RUNGE--KUTTA公式，等步长计算NSTEP步 SUBROUTINE RKDUMB(VSTART，NVAR，X1，X2，NSTEP，DERIVS) PARAMETER (NMAX=10) 函数最大个数 COMMON /PATH/ XX(200)，Y(10,200) DIMENSION VSTART(NVAR)， V(NMAX)， DV(NMAX) 置初值 DO 13 K=1，NSTEP CALL DERIVS(X，V，DV) CALL RK4(V，DV，,NVAR，X，H，V，DERIVS) X=X+H XX(K+1)=X 存贮中间步骤 DO 12 I=1，NVAR Y(I，K+1)=V(I) CONTINUE 结束返回 12 13 浙江大学研究生学位课程 《实用数值计算方法》
带自适应步长控制的四阶Runge--Kutta 驱动程序 6.2.4 带自适应步长控制的四阶Runge--Kutta 驱动程序 SUBROUTINE ODEINT(YSTART，NVAR，X1，X2，EPS，H1 HMIN，NOK，NBAD，DERIVS,，RKQC) PARAMETER (MAXSTEP=10000,NMAX=10,TWO=2.0, ZERO=0.0,TINY=1.0E-30) COMMON /PATH/ KMAX，KOUNT，DXSAV，XP(200)，YP(10,200) DIMENSION YSTART(NVAR)，YSCAL(NMAX)，Y(NMAX)，DYDX(NMAX) X=X1 H=SIGH(H1，X2-X1) NOK=0 NBAD=0 KOUNT=0 DO 11 I=1，NVAR Y(I)=YSTART(I) XSAV=X-DXSAV TWO ----保证第一步存贮 DO 16 NSTEP=1, MAXSTEP ----最多取MAXSTEP步数 CALL DERIVS(X，Y，DYDX) DO 12 I=1，NVAR 改变比例因子控制准确性 YSCAL(I)=ABS(Y(I))+ABS(H DYDX(I)) + TINY IF(KMAX.GT.O.AND.ABS(X-XSAV).GT.ABS(DXSAV) KOUNT. LT. KMAX-1) THEN KOUNT=KOUNT+1 XP(KOUNT)=X DO 13 I=1，NVAR YP(I，KOUNT)=Y(I) XSAV=X ENDIF 置初值 11 12 存贮中间结果 13 浙江大学研究生学位课程 《实用数值计算方法》
6.2.4 IF ((X+H-X2) (X+H-X1) .GT. ZERO) H=X2-X1 越界处理 CALL RKQC(Y，DYDX，NVAR，X,H,EPS,YSCAL,HDID,HNEXT,DERIVS) IF (HDID. EQ. H) THEN NOK=NOK+1 好步长计算 ELSE NBAD=NBAD+1 坏步长计算 ENDIF IF ((X-X2) (X2-X1). GE. ZERO) THEN 进行完毕 DO 14 I=1，NVAR YSTART(I)=Y(I) IF (KMAX. NE. O) THEN 保留最终步结果 KOUNT=KOUNT+1 XP(KOUNT)=X DO 15 I=1，NVAR YP(I，KOUNT)=Y(I) RERURN IF (ABS(HNEXT). LT. HMIN) PAUSE '步长太小' H=HNEXT CONTINUE PAUSE '步数太多，越界' RETURN END 14 15 16 浙江大学研究生学位课程 《实用数值计算方法》
6.2.4 浙江大学研究生学位课程 《实用数值计算方法》
（Stiff Sets of Equation） 6.2.5 方程组与刚性问题 （Stiff Sets of Equation） 1. 考虑2阶常微分方程组的初值问题： 解此方程组的Euler方法 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 可能是复数 浙江大学研究生学位课程 《实用数值计算方法》
绝对收敛。 Adams内插法的绝对稳定区域： 图 6.11 Adams内插法的绝对稳定区域 6.2.5 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 图6.12 Adams外推法的绝对稳定区域 图 6.13 R-K法的绝对稳定区域 浙江大学研究生学位课程 《实用数值计算方法》
2. 刚性方程组 （Stiff Equations） 6.2.5 2. 刚性方程组 （Stiff Equations） 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 定义： 刚性方程： 刚性比： 一般坏条件问题 严重坏条件问题 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 图 6.14 A-稳定区域 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 图 6.15 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 浙江大学研究生学位课程 《实用数值计算方法》
6.2.5 浙江大学研究生学位课程 《实用数值计算方法》
6.2.4 浙江大学研究生学位课程 《实用数值计算方法》
Euler法 一阶 Euler 预估--校正法 二阶 中点法二阶 (Runge--Kutta) 四级四阶 RungeKutta 6.2.4 Euler法 一阶 Euler 预估--校正法 二阶 中点法二阶 (Runge--Kutta) 四级四阶 RungeKutta 浙江大学研究生学位课程 《实用数值计算方法》
6.3 边值问题解法 浙江大学研究生学位课程 《实用数值计算方法》
6.3.1 试射法（Shooting） 图 6.16 试射法求解示意 浙江大学研究生学位课程 《实用数值计算方法》
6.3.1 浙江大学研究生学位课程 《实用数值计算方法》
6.3.1 浙江大学研究生学位课程 《实用数值计算方法》
6.3.2 差分方法（Difference Methods） (Relaxation Methods) Initial guess 1st iteration 2nd iteration True solution 图 6.17 差分方法示意 浙江大学研究生学位课程 《实用数值计算方法》
6.3.2 浙江大学研究生学位课程 《实用数值计算方法》
6.3.2 浙江大学研究生学位课程 《实用数值计算方法》
6.3.2 图 6.18 修正值矩阵 浙江大学研究生学位课程 《实用数值计算方法》
利用矩阵的特殊结构使总运算次数达到极小， 并使矩阵系数的存贮达到极小 6.3.2 块对角线性代数方程组的快速求解 利用矩阵的特殊结构使总运算次数达到极小， 并使矩阵系数的存贮达到极小 图 6.19 Gauss 消去法的目标结构 浙江大学研究生学位课程 《实用数值计算方法》
a) b) D--有待于对角化 A--将要改动 S--需要存贮 Z--将要变为0 图 6.20 特殊矩阵结构 1 6.3.2 图 6.20 特殊矩阵结构 1 浙江大学研究生学位课程 《实用数值计算方法》
a) b) 包括四步运算 1. 使用前一步的结果，简化部分元素为零 2. 消去剩下子块正方结构为单位阵 3. 存储以后要使用的非零系数 6.3.2 a) b) 图 6.22 特殊矩阵结构 2 包括四步运算 1. 使用前一步的结果，简化部分元素为零 2. 消去剩下子块正方结构为单位阵 3. 存储以后要使用的非零系数 4. 回代 浙江大学研究生学位课程 《实用数值计算方法》
6.3.3 样条函数在两点边值问题中的应用 考察 二阶线性微分方程 的边值问题 浙江大学研究生学位课程 《实用数值计算方法》
6.3.3 浙江大学研究生学位课程 《实用数值计算方法》
6.3.3 浙江大学研究生学位课程 《实用数值计算方法》
6.3.3 浙江大学研究生学位课程 《实用数值计算方法》
6.3.3 浙江大学研究生学位课程 《实用数值计算方法》
6.3.3 浙江大学研究生学位课程 《实用数值计算方法》
1914年--1923年捕获的鱼中，食用鱼所占的比例(%) 6.3.3 相平面上的轨迹为 图 6.23 相平面上的轨迹图 1914年--1923年捕获的鱼中，食用鱼所占的比例(%) 浙江大学研究生学位课程 《实用数值计算方法》
6.3.3 浙江大学研究生学位课程 《实用数值计算方法》
总 结 公式简单，但精度不好 编程简单，使用方便 计算量少，但稳定性差 计算量大，但稳定性好 6.3.3 浙江大学研究生学位课程 总 结 公式简单，但精度不好 编程简单，使用方便 计算量少，但稳定性差 计算量大，但稳定性好 浙江大学研究生学位课程 《实用数值计算方法》
第六章 习题 浙江大学研究生学位课程 《实用数值计算方法》
浙江大学研究生学位课程 《实用数值计算方法》