【教学课件】第六章常微分方程及方程组的解法.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,数值解及其重要性,浙江大学研究生学位课程,9,ODE数值解的基本思想和方法特点 基本思想有两点 1.离散化 用Taylor级数，数值积分和差商 逼近导数等手段，把ODE转化为离散 的代数方程(称差分方程)。2.递推化 在具有唯一解的条件下，通过 步进法逐步计算出解在一系列离散 点上的值。从而得到原ODE的数值 近似解。,浙江大学研究生学位课程,10,6.2 初值问题解法 我们讨论一阶ODE，而高阶可 能化为一阶ODEs。一阶初值问题 可以一般地写成：,欧拉（Euler）方法 Euler方法是求解(6-8)最简单方法，但精度差，故不实用。然而对理论分 析很有用。,浙江大学研究生学位课程,11,6.2.1.1 方法原理及推导 设初值问题(6-8)满足：,浙江大学研究生学位课程,12,6.2.1,图 6.1 常微分方程初值问题的数值解,浙江大学研究生学位课程,13,浙江大学研究生学位课程,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公式求解,浙江大学研究生学位课程,19,图 6.4 运动轨迹,浙江大学研究生学位课程,20,图 6.5 Euler公式的误差,6.2.1.2 Euler方法的误差估计 一般其它方法的误差估计也类似。这里误差是指截断误差(算法理论误 差)而不是舍入误差。后者由计算机字 长等决定，属于稳定性问题。i)几何分析,6.2.1,浙江大学研究生学位课程,21,6.2.1,局部截断,误差，,浙江大学研究生学位课程,22,整体截断误差：设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,浙江大学研究生学位课程,30,6.2.1,浙江大学研究生学位课程,31,浙江大学研究生学位课程,32,表 6.2 予估校正求解结果对比,浙江大学研究生学位课程,33,表 6.3 Euler法与外推结果的比较,浙江大学研究生学位课程,34,线性多步法,浙江大学研究生学位课程,35,(),(),浙江大学研究生学位课程,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,预估-校正方法（Predictor-Corrector Method）,浙江大学研究生学位课程,53,6.2.2,浙江大学研究生学位课程,54,注意：一步校正的计算量 预估计算量。所以要适当选取h才能发挥PC的优点。设绝对稳定区域：达到精度的校正次数为N，则h的选取，应满足：否则可用N步显式算法稳定达到目的。,h,浙江大学研究生学位课程,55,Adams四步四阶预估校正算法,浙江大学研究生学位课程,56,6.2.3 Runge-Kutta 方法,浙江大学研究生学位课程,57,6.2.3,浙江大学研究生学位课程,58,6.2.3,浙江大学研究生学位课程,59,6.2.3,浙江大学研究生学位课程,60,6.2.3,六个未知数，二个自由，故可取,故：,浙江大学研究生学位课程,61,N=4:四级四阶R-K方法,-最常用的古典Runge-Kutta方法。增加计算函数值的次数（级）与提高精 度（阶）的关系见下表：,表 6.7 Runge-Kutta方法中级与阶的关系,浙江大学研究生学位课程,62,6.2.3,表 6.8 各种解法在例题中的结果比较,浙江大学研究生学位课程,63,（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-4EXTERNAL DERIVSDIMENSION 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=HTRYHH=0.5HCALL RK4(YSAV,DYSAV,N,XSAV,HH,YTEMP,DERIVS)X=XSAV+HHCALL DERIVS(X,YTEMP,DYDX)CALL RK4(YTMP,DYDX,N,X,HH,Y,DERIVS)X=XSAV+HIF(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 RKQCX1=1.0X2=10.0VS(1)=BESSO(X1)VS(2)=BESSI(X1)VS(3)=BESSJ(2,X1)VS(4)=BESSJ(3,X1)EPS=1.0E-4HI=1.0HMIN=0.0KMAX=100DXSAV=(X1-X2)/20.0CALL ODEINT(VS，NVAR，X1，X2，EPS，HI，HMIN，NOK NBAD，NOK，NBAD，DERIVS，RKQC)WRITE(，)ENDSUBROUTINE 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)RETURNEND,6.2.4,浙江大学研究生学位课程,74,IF(ERRMAX.GT.ONE)THEN 不满足要求 H=SAFETY H(ERRMAX PSHRINK)估计下次步长（缩小）GOTO 1ELSE 满足要求 HDID=H IF(ERRMAX.GT.ERRCON)THEN HNEXT=SAFETY H(ERRMAX PGROW)估计可放大的步长 ELSE HNEXT=4.H 最多可放大四倍 ENDIFENDIFDO 13 I=1，N 完成五阶截断误差 Y(I)=Y(I)+YTEMP(I)FCORCONTINUERETURNEND四阶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.5H6=H/6.0XH=X+HHDO 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)CONTINUERETURNEND,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)CONTINUECONTINUE结束返回,1213,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=X1H=SIGH(H1，X2-X1)NOK=0NBAD=0KOUNT=0DO 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 坏步长计算ENDIFIF(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 RERURNENDIFIF(ABS(HNEXT).LT.HMIN)PAUSE 步长太小H=HNEXTCONTINUEPAUSE 步数太多，越界RETURNEND,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.12 Adams外推法的绝对稳定区域,图 6.13 R-K法的绝对稳定区域,浙江大学研究生学位课程,86,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,浙江大学研究生学位课程,92,6.2.5,浙江大学研究生学位课程,93,浙江大学研究生学位课程,94,6.2.4,浙江大学研究生学位课程,95,Euler法一阶,中点法二阶(Runge-Kutta),四级四阶RungeKutta,Euler 预估-校正法二阶,6.2.4,浙江大学研究生学位课程,96,6.3 边值问题解法,浙江大学研究生学位课程,97,试射法（Shooting）,图 6.16 试射法求解示意,浙江大学研究生学位课程,98,浙江大学研究生学位课程,99,浙江大学研究生学位课程,100,6.3.2 差分方法（Difference Methods）(Relaxation Methods),True solution,2nd iteration,1st iteration,Initial guess,图 6.17 差分方法示意,浙江大学研究生学位课程,101,6.3.2,浙江大学研究生学位课程,102,浙江大学研究生学位课程,103,图 6.18 修正值矩阵,浙江大学研究生学位课程,104,块对角线性代数方程组的快速求解 利用矩阵的特殊结构使总运算次数达到极小，并使矩阵系数的存贮达到极小,图 6.19 Gauss 消去法的目标结构,浙江大学研究生学位课程,105,D-有待于对角化A-将要改动S-需要存贮Z-将要变为0,a)b)a)b),图 6.20 特殊矩阵结构 1,浙江大学研究生学位课程,106,a)b),包括四步运算1.使用前一步的结果，简化部分元素为零2.消去剩下子块正方结构为单位阵3.存储以后要使用的非零系数4.回代,图 6.22 特殊矩阵结构 2,浙江大学研究生学位课程,107,样条函数在两点边值问题中的应用 考察 二阶线性微分方程 的边值问题,浙江大学研究生学位课程,108,浙江大学研究生学位课程,109,6.3.3,浙江大学研究生学位课程,110,6.3.3,浙江大学研究生学位课程,111,浙江大学研究生学位课程,112,浙江大学研究生学位课程,113,相平面上的轨迹为,1914年-1923年捕获的鱼中，食用鱼所占的比例(%),图 6.23 相平面上的轨迹图,浙江大学研究生学位课程,114,浙江大学研究生学位课程,115,公式简单，但精度不好,编程简单，使用方便,计算量少，但稳定性差,计算量大，但稳定性好,总 结,浙江大学研究生学位课程,116,第六章 习题,浙江大学研究生学位课程,117,