精品]浙江大学研究生学位课程.ppt

浙江大学研究生学位课程,实用数值计算方法,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 基本概念 描述自由落体的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，如果能写成：,则称该方程为显式的。方程（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级数，数值积分和差商 逼近导数等手段，把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,浙江大学研究生学位课程,实用数值计算方法,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.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）的解。 则 整体截断误差定义为 它是局部截断误差的积累。 定理：若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,浙江大学研究生学位课程,实用数值计算方法,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步外插法,浙江大学研究生学位课程,实用数值计算方法,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,浙江大学研究生学位课程,实用数值计算方法,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 预估-校正方法 （Predictor-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,浙江大学研究生学位课程,实用数值计算方法,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.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.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.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 步长无意义 CALL RK4 (YSAV, DYSAV, N, XSAV, HH, YTEMP, DERIVS) 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,6.2.4,浙江大学研究生学位课程,实用数值计算方法,73,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 END,6.2.4,浙江大学研究生学位课程,实用数值计算方法,74,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 ENDIF DO 13 I=1，N 完成五阶截断误差 Y(I)=Y(I) + YTEMP(I) FCOR CONTINUE RETURN END 四阶Runge-Kutta 步进程序（续）,13,6.2.4,浙江大学研究生学位课程,实用数值计算方法,75,计算结果 NOK=30 NBAD=0 KOUNT=15,6.2.4,浙江大学研究生学位课程,实用数值计算方法,76,四阶 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,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 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 CONTINUE 结束返回,12 13,6.2.4,浙江大学研究生学位课程,实用数值计算方法,78,带自适应步长控制的四阶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,存贮中间结果,13,12,11,置初值,6.2.4,浙江大学研究生学位课程,实用数值计算方法,79,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) ENDIF RERURN ENDIF 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,可能是复数,浙江大学研究生学位课程,实用数值计算方法,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,定义：,刚性方程：,一般坏条件问题,严重坏条件问题,刚性比：,浙江大学研究生学位课程,实用数值计算方法,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.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 iteratio