1、浙江大学研究生学位课程,实用数值计算方法,1,第六章 常微分方程及 方程组的解法,6.1 常微分方程及其求解概述 6.2 初值问题解法 6.3 边值问题解法,浙江大学研究生学位课程,实用数值计算方法,2,6.1 常微分方程及其求解概述 6.2 初值问题解法1.Euler方法2.线性多步法3.Runge-Kutta法4.方程组及刚性问题的Gear方法 6.3 边值问题解法1.Shooting(试射法)2.差分法,浙江大学研究生学位课程,实用数值计算方法,3,6.1 常微分方程及求解概述 （Ordinary Differential Equations, ODE),6.1.1 基本概念描述自由落体
2、的ODE:,浙江大学研究生学位课程,实用数值计算方法,4, 只有一个自变量的微分方程为 ODE,否则称为偏微分方程 PDE。 方程中未知函数导数的最高阶数称为方程的阶。 （6-4）是二阶的 方程中关于未知函数及其各阶导数均是一次的，则称为线性微分方程。,和（6-1）都是线性二阶ODE。 (6-2),(6-3)是(6-1)的初始条件。亦称定解条件。(6-1)(6-2)叫做初值问题。,6.1.1,浙江大学研究生学位课程,实用数值计算方法,5,6.1.1,（6-1），（6-3）叫做边值问题。 在没有给定解条件时。方程一般有一族解曲线 y(x,c) 。如：, 对任意的n阶ODE，如果能写成：,则称该方
3、程为显式的。方程（6-4）是显式的。而下面方程是隐式的。,浙江大学研究生学位课程,实用数值计算方法,6, 对于高阶显式方程。通过定义n-1个新变量，可以写成n维一阶方程组。即令：,6.1.1,浙江大学研究生学位课程,实用数值计算方法,7, 在讨论初值问题时，我们从一阶方程开始：,然后毫不费力地套用来解方程组。 当 f(x,y)与y无关时，f(x,y)=g(x),6.1.1,浙江大学研究生学位课程,实用数值计算方法,8,6.1.2 数值解及其重要性,浙江大学研究生学位课程,实用数值计算方法,9,6.1.3 ODE数值解的基本思想和方法特点基本思想有两点1. 离散化用Taylor级数，数值积分和差
4、商逼近导数等手段，把ODE转化为离散的代数方程(称差分方程)。2. 递推化在具有唯一解的条件下，通过步进法逐步计算出解在一系列离散点上的值。从而得到原ODE的数值近似解。,浙江大学研究生学位课程,实用数值计算方法,10,6.2 初值问题解法我们讨论一阶ODE，而高阶可能化为一阶ODEs。一阶初值问题可以一般地写成：,6.2.1 欧拉（Euler）方法Euler方法是求解(6-8)最简单方法，但精度差，故不实用。然而对理论分 析很有用。,浙江大学研究生学位课程,实用数值计算方法,11,6.2.1.1 方法原理及推导设初值问题(6-8)满足：,6.2.1,浙江大学研究生学位课程,实用数值计算方法,
5、12,6.2.1,图 6.1 常微分方程初值问题的数值解,浙江大学研究生学位课程,实用数值计算方法,13,6.2.1,浙江大学研究生学位课程,实用数值计算方法,14,欧拉方法的几何意义：,h步长,6.2.1,图 6.2 Euler方法的几何意义,浙江大学研究生学位课程,实用数值计算方法,15,6.2.1,浙江大学研究生学位课程,实用数值计算方法,16,6.2.1,浙江大学研究生学位课程,实用数值计算方法,17,6.2.1,多 步 法,单步法,自动起步,显式,隐式,半隐式,图 6.3 ODE求解方法的类型,浙江大学研究生学位课程,实用数值计算方法,18,表 6.1 自由落体运动方程的Euler公
6、式求解,6.2.1,浙江大学研究生学位课程,实用数值计算方法,19,图 6.4 运动轨迹,6.2.1,浙江大学研究生学位课程,实用数值计算方法,20,图 6.5 Euler公式的误差,6.2.1.2 Euler方法的误差估计一般其它方法的误差估计也类似。这里误差是指截断误差(算法理论误差)而不是舍入误差。后者由计算机字长等决定，属于稳定性问题。i) 几何分析,6.2.1,浙江大学研究生学位课程,实用数值计算方法,21,6.2.1,局部截断,误差，,浙江大学研究生学位课程,实用数值计算方法,22,6.2.1 整体截断误差：设ym是Euler公式（6-9）精确解，而 y(x) 是初值问题（6-8）
7、的解。则 整体截断误差定义为它是局部截断误差的积累。定理：若f(x,y)关于y满足Lipschitz条件。则有估计式：,浙江大学研究生学位课程,实用数值计算方法,23,6.2.1,浙江大学研究生学位课程,实用数值计算方法,24,6.2.1,浙江大学研究生学位课程,实用数值计算方法,25,6.2.1,浙江大学研究生学位课程,实用数值计算方法,26,6.2.1,注意：,稳定性：,浙江大学研究生学位课程,实用数值计算方法,27,6.2.1,浙江大学研究生学位课程,实用数值计算方法,28,6.2.1,浙江大学研究生学位课程,实用数值计算方法,29,6.2.1,浙江大学研究生学位课程,实用数值计算方法,
8、30,6.2.1,浙江大学研究生学位课程,实用数值计算方法,31,6.2.1,浙江大学研究生学位课程,实用数值计算方法,32,6.2.1,表 6.2 予估校正求解结果对比,浙江大学研究生学位课程,实用数值计算方法,33,6.2.1,表 6.3 Euler法与外推结果的比较,浙江大学研究生学位课程,实用数值计算方法,34,6.2.2 线性多步法,浙江大学研究生学位课程,实用数值计算方法,35,6.2.2,(),(),浙江大学研究生学位课程,实用数值计算方法,36,6.2.2,Adams 外插法 (k=2) 3阶3步 显式,表 6.4 外插系数bki值,图 6.6 3阶3步外插法,浙江大学研究生学
9、位课程,实用数值计算方法,37,6.2.2,浙江大学研究生学位课程,实用数值计算方法,38,6.2.2,Adams 外插法 (k=2) 4阶3步,图 6.7 4步3阶Adams内插公式,浙江大学研究生学位课程,实用数值计算方法,39,6.2.2,浙江大学研究生学位课程,实用数值计算方法,40,6.2.2,浙江大学研究生学位课程,实用数值计算方法,41,6.2.2,图 6.8 一般化插值形式,浙江大学研究生学位课程,实用数值计算方法,42,6.2.2,浙江大学研究生学位课程,实用数值计算方法,43,6.2.2,浙江大学研究生学位课程,实用数值计算方法,44,6.2.2,浙江大学研究生学位课程,实
10、用数值计算方法,45,6.2.2,浙江大学研究生学位课程,实用数值计算方法,46,6.2.2,浙江大学研究生学位课程,实用数值计算方法,47,6.2.2,图 6.9,浙江大学研究生学位课程,实用数值计算方法,48,6.2.2,线性多步法的绝对稳定性：,浙江大学研究生学位课程,实用数值计算方法,49,6.2.2,定义：,绝对稳定。,绝对稳定区域。,浙江大学研究生学位课程,实用数值计算方法,50,6.2.2,Milne,浙江大学研究生学位课程,实用数值计算方法,51,6.2.2,表 6.6 计算结果,浙江大学研究生学位课程,实用数值计算方法,52,6.2.2,6.2.2.4 预估-校正方法（Pre
11、dictor-Corrector Method）,浙江大学研究生学位课程,实用数值计算方法,53,6.2.2,浙江大学研究生学位课程,实用数值计算方法,54,6.2.2 注意：一步校正的计算量 预估计算量。所以要适当选取h才能发挥PC的优点。设绝对稳定区域：达到精度的校正次数为N，则h的选取，应满足：否则可用N步显式算法稳定达到目的。,h,浙江大学研究生学位课程,实用数值计算方法,55,6.2.2 Adams四步四阶预估校正算法,浙江大学研究生学位课程,实用数值计算方法,56,6.2.3 Runge-Kutta 方法,浙江大学研究生学位课程,实用数值计算方法,57,6.2.3,浙江大学研究生学
12、位课程,实用数值计算方法,58,6.2.3,浙江大学研究生学位课程,实用数值计算方法,59,6.2.3,浙江大学研究生学位课程,实用数值计算方法,60,6.2.3,六个未知数，二个自由，故可取,故：,浙江大学研究生学位课程,实用数值计算方法,61,6.2.3 N=4: 四级四阶R-K方法,-最常用的古典Runge-Kutta方法。增加计算函数值的次数（级）与提高精度（阶）的关系见下表：,表 6.7 Runge-Kutta方法中级与阶的关系,浙江大学研究生学位课程,实用数值计算方法,62,6.2.3,表 6.8 各种解法在例题中的结果比较,浙江大学研究生学位课程,实用数值计算方法,63,6.2.
13、3 （2）. 单步法，自动起步（3）. 易改为变步长（4）. 绝对稳定区域较同阶线性多步法大（5）. 计算工作量较大，有时大于隐式方法（6）. 估计误差不易绝对稳定性讨论,浙江大学研究生学位课程,实用数值计算方法,64,6.2.3,表 6.9 各级R-K方法的绝对稳定区域,浙江大学研究生学位课程,实用数值计算方法,65,6.2.3,表 6.10 不同步长对精度的影响,浙江大学研究生学位课程,实用数值计算方法,66,6.2.3,浙江大学研究生学位课程,实用数值计算方法,67,6.2.3,浙江大学研究生学位课程,实用数值计算方法,68,P阶 图 6.10 变步长Runge-Kutta方法框图,6.
14、2.3,浙江大学研究生学位课程,实用数值计算方法,69,6.2.3,浙江大学研究生学位课程,实用数值计算方法,70,6.2.4 出发值的计算,浙江大学研究生学位课程,实用数值计算方法,71,6.2.4,浙江大学研究生学位课程,实用数值计算方法,72,质量控制 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
15、.E-4 EXTERNAL DERIVS DIMENSION Y(N), DYSX(N), YSCAL(N), YTEMP(NMAX), YSAV(NMAX), DYSAV(NMAX) XSAV=X 保留初值 DO 11 I=1，NYSAV(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+
16、H IF (X. EQ. XSAV) PAUS E 步长无意义 CALL RK4 (YSAV, DYSAV, N, XSAV, HH, YTEMP, DERIVS) ERRMAX=0. DO 12 I=1，NYTEMP(I)=Y(I)-YTEMP(I)ERRMAX = MAX(ERRMAX, ABS(YTEMP(I) / YSCAL(I) ERRMAX = ERRMAX / EPS,误差计算,一个单步计算,两个半步长计算,置初始值,11,1,12,6.2.4,浙江大学研究生学位课程,实用数值计算方法,73,PARAMETER (NVAR=4) DIMENSION VS(NVAR) COMMON
17、 /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，NOKNBAD， NOK，NBAD，DERIVS，RKQC) WRITE ( ， ) END SUBROUTINE DE
18、RIVS (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 END,6.2.4,浙江大学研究生学位课程,实用数值计算方法,74,IF (ERRMAX.GT.ONE) THEN 不满足要求H = SAFETY H (ERRMAX PSHRINK) 估计下次步长（缩小）GOTO 1 ELSE 满足要求HDID = HIF (ERRMAX. GT. ERRCON) THE
19、NHNEXT = SAFETY H (ERRMAX PGROW) 估计可放大的步长ELSEHNEXT =4. H 最多可放大四倍ENDIF ENDIF DO 13 I=1，N 完成五阶截断误差Y(I)=Y(I) + YTEMP(I) FCOR CONTINUE RETURN END四阶Runge-Kutta 步进程序（续）,13,6.2.4,浙江大学研究生学位课程,实用数值计算方法,75,计算结果NOK=30NBAD=0KOUNT=15,6.2.4,浙江大学研究生学位课程,实用数值计算方法,76,四阶 Runge-Kutta 方法： SUBROUTINE RK4(Y，DYDX，N，X，H，YO
20、UT，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，NYT(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 DER
21、IVS(X+H，YT，DYT) DO 14 I=1，NYOUT(I)=DYDX(I)+DYT(I)+2. DYM(I)YOUT(I)=Y(I) + H6 YOUT(I) CONTINUE RETURN END,14,13,12,11,6.2.4,浙江大学研究生学位课程,实用数值计算方法,77,用四阶 RUNGE-KUTTA公式，等步长计算NSTEP步 SUBROUTINE RKDUMB(VSTART，NVAR，X1，X2，NSTEP，DERIVS) PARAMETER (NMAX=10) 函数最大个数 COMMON /PATH/ XX(200)，Y(10,200) DIMENSION VSTA
22、RT(NVAR)， V(NMAX)， DV(NMAX)置初值DO 13 K=1，NSTEPCALL DERIVS(X，V，DV)CALL RK4(V，DV，,NVAR，X，H，V，DERIVS)X=X+HXX(K+1)=X 存贮中间步骤DO 12 I=1，NVARY(I，K+1)=V(I)CONTINUE CONTINUE结束返回,12 13,6.2.4,浙江大学研究生学位课程,实用数值计算方法,78,带自适应步长控制的四阶Runge-Kutta 驱动程序 SUBROUTINE ODEINT(YSTART，NVAR，X1，X2，EPS，H1HMIN，NOK，NBAD，DERIVS,，RKQC)
23、 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，NVARY(I)=YSTART(I) XSAV=X-DXSAV TWO -保证第一步存贮 DO 16 NSTEP=1, MAXSTEP -最多取MAXSTE
24、P步数CALL DERIVS(X，Y，DYDX)DO 12 I=1，NVAR 改变比例因子控制准确性YSCAL(I)=ABS(Y(I)+ABS(H DYDX(I) + TINYIF(KMAX.GT.O.AND.ABS(X-XSAV).GT.ABS(DXSAV)KOUNT. LT. KMAX-1) THENKOUNT=KOUNT+1XP(KOUNT)=XDO 13 I=1，NVARYP(I，KOUNT)=Y(I)XSAV=XENDIF,存贮中间结果,13,12,11,置初值,6.2.4,浙江大学研究生学位课程,实用数值计算方法,79,IF (X+H-X2) (X+H-X1) .GT. ZERO)
25、 H=X2-X1 越界处理 CALL RKQC(Y，DYDX，NVAR，X,H,EPS,YSCAL,HDID,HNEXT,DERIVS) IF (HDID. EQ. H) THENNOK=NOK+1 好步长计算 ELSENBAD=NBAD+1 坏步长计算 ENDIF IF (X-X2) (X2-X1). GE. ZERO) THEN 进行完毕DO 14 I=1，NVARYSTART(I)=Y(I)IF (KMAX. NE. O) THEN 保留最终步结果KOUNT=KOUNT+1XP(KOUNT)=XDO 15 I=1，NVARYP(I，KOUNT)=Y(I)ENDIFRERURN ENDIF
26、 IF (ABS(HNEXT). LT. HMIN) PAUSE 步长太小 H=HNEXT CONTINUE PAUSE 步数太多，越界 RETURN END,16,15,14,6.2.4,浙江大学研究生学位课程,实用数值计算方法,80,6.2.4,浙江大学研究生学位课程,实用数值计算方法,81,6.2.5 方程组与刚性问题（Stiff Sets of Equation）1. 考虑2阶常微分方程组的初值问题：,解此方程组的Euler方法,浙江大学研究生学位课程,实用数值计算方法,82,6.2.5,浙江大学研究生学位课程,实用数值计算方法,83,6.2.5,可能是复数,浙江大学研究生学位课程,实
27、用数值计算方法,84,6.2.5,绝对收敛。,Adams内插法的绝对稳定区域：,图 6.11 Adams内插法的绝对稳定区域,浙江大学研究生学位课程,实用数值计算方法,85,6.2.5,图6.12 Adams外推法的绝对稳定区域,图 6.13 R-K法的绝对稳定区域,浙江大学研究生学位课程,实用数值计算方法,86,6.2.5 2. 刚性方程组 （Stiff Equations）,浙江大学研究生学位课程,实用数值计算方法,87,6.2.5,浙江大学研究生学位课程,实用数值计算方法,88,6.2.5,定义：,刚性方程：,一般坏条件问题,严重坏条件问题,刚性比：,浙江大学研究生学位课程,实用数值计算
28、方法,89,6.2.5,图 6.14 A-稳定区域,浙江大学研究生学位课程,实用数值计算方法,90,6.2.5,图 6.15,浙江大学研究生学位课程,实用数值计算方法,91,6.2.5,浙江大学研究生学位课程,实用数值计算方法,92,6.2.5,浙江大学研究生学位课程,实用数值计算方法,93,6.2.5,浙江大学研究生学位课程,实用数值计算方法,94,6.2.4,浙江大学研究生学位课程,实用数值计算方法,95,Euler法 一阶,中点法二阶 (Runge-Kutta),四级四阶 RungeKutta,Euler 预估-校正法 二阶,6.2.4,浙江大学研究生学位课程,实用数值计算方法,96,6
29、.3 边值问题解法,浙江大学研究生学位课程,实用数值计算方法,97,6.3.1 试射法（Shooting）,图 6.16 试射法求解示意,浙江大学研究生学位课程,实用数值计算方法,98,6.3.1,浙江大学研究生学位课程,实用数值计算方法,99,6.3.1,浙江大学研究生学位课程,实用数值计算方法,100,6.3.2 差分方法（Difference Methods）(Relaxation Methods),True solution,2nd iteration,1st iteration,Initial guess,图 6.17 差分方法示意,浙江大学研究生学位课程,实用数值计算方法,101,
30、6.3.2,浙江大学研究生学位课程,实用数值计算方法,102,6.3.2,浙江大学研究生学位课程,实用数值计算方法,103,6.3.2,图 6.18 修正值矩阵,浙江大学研究生学位课程,实用数值计算方法,104,块对角线性代数方程组的快速求解利用矩阵的特殊结构使总运算次数达到极小，并使矩阵系数的存贮达到极小,图 6.19 Gauss 消去法的目标结构,6.3.2,浙江大学研究生学位课程,实用数值计算方法,105,D-有待于对角化 A-将要改动 S-需要存贮 Z-将要变为0,a)b)a)b),6.3.2,图 6.20 特殊矩阵结构 1,浙江大学研究生学位课程,实用数值计算方法,106,a)b),
31、包括四步运算 1. 使用前一步的结果，简化部分元素为零 2. 消去剩下子块正方结构为单位阵 3. 存储以后要使用的非零系数 4. 回代,6.3.2,图 6.22 特殊矩阵结构 2,浙江大学研究生学位课程,实用数值计算方法,107,6.3.3 样条函数在两点边值问题中的应用考察 二阶线性微分方程 的边值问题,浙江大学研究生学位课程,实用数值计算方法,108,6.3.3,浙江大学研究生学位课程,实用数值计算方法,109,6.3.3,浙江大学研究生学位课程,实用数值计算方法,110,6.3.3,浙江大学研究生学位课程,实用数值计算方法,111,6.3.3,浙江大学研究生学位课程,实用数值计算方法,112,6.3.3,浙江大学研究生学位课程,实用数值计算方法,113,相平面上的轨迹为,1914年-1923年捕获的鱼中，食用鱼所占的比例(%),图 6.23 相平面上的轨迹图,6.3.3,浙江大学研究生学位课程,实用数值计算方法,114,6.3.3,浙江大学研究生学位课程,实用数值计算方法,115,公式简单，但精度不好,编程简单，使用方便,计算量少，但稳定性差,计算量大，但稳定性好,6.3.3,总 结,浙江大学研究生学位课程,实用数值计算方法,116,第六章 习题,浙江大学研究生学位课程,实用数值计算方法,117,
