六常微分方程及方程组的解法.pptx

浙江大学研究生学位课程16.1 常微分方程及其求解概述常微分方程及其求解概述6.2 初值问题解法初值问题解法 1.Euler方法方法 2.线性多步法线性多步法 3.Runge-Kutta法法 4.方程组及刚性问题的方程组及刚性问题的Gear方法方法6.3 边值问题解法边值问题解法 1.Shooting(试射法试射法)2.差分法差分法 浙江大学研究生学位课程26.1 常微分方程及求解概述常微分方程及求解概述（Ordinary Differential Equations,ODE)6.1.1 基本概念基本概念 描述自由落体的描述自由落体的ODE:浙江大学研究生学位课程3 只有一个自变量的微分方程为只有一个自变量的微分方程为 ODE,否则称为偏微分方程否则称为偏微分方程 PDE。方程中未知函数导数的最高阶数称为方程中未知函数导数的最高阶数称为 方程的阶。方程的阶。（6-4）是二阶的）是二阶的 方程中关于未知函数及其各阶导数方程中关于未知函数及其各阶导数 均是一次的，则称为线性微分方程。均是一次的，则称为线性微分方程。和（和（6-1）都是线性二阶）都是线性二阶ODE。(6-2),(6-3)是是(6-1)的的初始条件初始条件。亦称定解。亦称定解 条件。条件。(6-1)(6-2)叫做叫做初值问题初值问题。6.1.1 浙江大学研究生学位课程46.1.1 （6-1），（6-3）叫做边值问题。在没有给定解条件时。方程一般在没有给定解条件时。方程一般 有一族解曲线有一族解曲线 y(x,c)。如：。如：对任意的对任意的n阶阶ODE，如果能写成：，如果能写成：则称该方程为则称该方程为显式显式的。方程（的。方程（6-4）是）是 显式的。而下面方程是显式的。而下面方程是隐式隐式的。的。浙江大学研究生学位课程5 对于高阶显式方程。通过定义对于高阶显式方程。通过定义n-1个个 新变量，可以写成新变量，可以写成n维一阶方程组。维一阶方程组。即令：即令：6.1.1 浙江大学研究生学位课程6 在讨论初值问题时，我们从一阶方在讨论初值问题时，我们从一阶方 程开始：程开始：然后毫不费力地套用来解方程组。然后毫不费力地套用来解方程组。当当 f(x,y)与与y无关时，无关时，f(x,y)=g(x)6.1.1 浙江大学研究生学位课程7 6.1.2 数值解及其重要性数值解及其重要性浙江大学研究生学位课程8 6.1.3 ODE数值解的基本思想和方法特点数值解的基本思想和方法特点 基本思想有两点基本思想有两点 1.离散化离散化 用Taylor级数，数值积分和差商 逼近导数等手段，把ODE转化为离散 的代数方程(称差分方程)。2.递推化递推化 在具有唯一解的条件下，通过 步进法逐步计算出解在一系列离散 点上的值。从而得到原ODE的数值 近似解。浙江大学研究生学位课程9 6.2 初值问题解法初值问题解法 我们讨论一阶ODE，而高阶可 能化为一阶ODEs。一阶初值问题 可以一般地写成：6.2.1 欧拉（欧拉（Euler）方法方法 Euler方法是求解(6-8)最简单简单方法，但精度差精度差，故不实用不实用。然而对理论分 析很有用。浙江大学研究生学位课程10 6.2.1.1 方法原理及推导方法原理及推导 设初值问题设初值问题(6-8)满足：满足：6.2.1浙江大学研究生学位课程116.2.1 图图 6.1 常微分方程初值问题的数值解常微分方程初值问题的数值解浙江大学研究生学位课程126.2.1浙江大学研究生学位课程13欧拉方法的几何意义：欧拉方法的几何意义：h步长6.2.1 图图 6.2 Euler方法的几何意义方法的几何意义浙江大学研究生学位课程146.2.1 浙江大学研究生学位课程156.2.1 浙江大学研究生学位课程166.2.1 多多步步法法单步法单步法自动起步自动起步显式显式隐式隐式半隐式半隐式图图 6.3 ODE求解方法的类型求解方法的类型浙江大学研究生学位课程17表表 6.1 自由落体运动方程的自由落体运动方程的Euler公式求解公式求解6.2.1浙江大学研究生学位课程18图图 6.4 运动轨迹运动轨迹6.2.1浙江大学研究生学位课程19图图 6.5 Euler公式的误差公式的误差 6.2.1.2 Euler方法的误差估计方法的误差估计 一般其它方法的误差估计也类似。这里误差是指截断误差截断误差(算法理论误 差)而不是舍入误差舍入误差。后者由计算机字 长等决定，属于稳定性问题。i)几何分析几何分析6.2.1 浙江大学研究生学位课程206.2.1 局部截断局部截断误差，误差，浙江大学研究生学位课程216.2.1 整体截断误差：整体截断误差：设ym是Euler公式（6-9）精确解，而 y(x)是初值问题（6-8）的解。则 整体截断误差定义为 它是局部截断误差的积累。定理：定理：若f(x,y)关于y满足Lipschitz条件。则有估计式：浙江大学研究生学位课程226.2.1 浙江大学研究生学位课程236.2.1 浙江大学研究生学位课程246.2.1 浙江大学研究生学位课程256.2.1 注意注意：稳定性稳定性：浙江大学研究生学位课程266.2.1 浙江大学研究生学位课程276.2.1 浙江大学研究生学位课程286.2.1浙江大学研究生学位课程296.2.1 浙江大学研究生学位课程306.2.1浙江大学研究生学位课程316.2.1表表 6.2 予估校正求解结果对比予估校正求解结果对比浙江大学研究生学位课程326.2.1表表 6.3 Euler法与外推结果的比较法与外推结果的比较浙江大学研究生学位课程33 6.2.2 线性多步法线性多步法 浙江大学研究生学位课程346.2.2()()浙江大学研究生学位课程356.2.2 Adams 外插法(k=2)3阶3步 显式表表 6.4 外插系数外插系数bki值值图图 6.6 3阶阶3步外插法步外插法浙江大学研究生学位课程366.2.2 浙江大学研究生学位课程376.2.2 Adams 外插法(k=2)4阶3步 图图 6.7 4步步3阶阶Adams内插公式内插公式浙江大学研究生学位课程386.2.2 浙江大学研究生学位课程396.2.2 浙江大学研究生学位课程406.2.2 图图 6.8 一般化插值形式一般化插值形式浙江大学研究生学位课程416.2.2 浙江大学研究生学位课程426.2.2 浙江大学研究生学位课程436.2.2 浙江大学研究生学位课程446.2.2 浙江大学研究生学位课程456.2.2 浙江大学研究生学位课程466.2.2 图图 6.9浙江大学研究生学位课程476.2.2 线性多步法的绝对稳定性：线性多步法的绝对稳定性：浙江大学研究生学位课程486.2.2 定义定义：绝对稳定绝对稳定。绝对稳定区域。绝对稳定区域。浙江大学研究生学位课程496.2.2 Milne浙江大学研究生学位课程506.2.2 表表 6.6 计算结果计算结果浙江大学研究生学位课程516.2.2 6.2.2.4 预估预估-校正方法校正方法 （Predictor-Corrector Method）浙江大学研究生学位课程526.2.2 浙江大学研究生学位课程536.2.2 注意：注意：一步校正的计算量 预估计算量。所以要适当选取h才能发挥PC的优点。设绝对稳定区域：设绝对稳定区域：达到精度的校正次数为达到精度的校正次数为N，则，则h的选取，的选取，应满足：应满足：否则可用否则可用N步显式算法稳定达到目的步显式算法稳定达到目的。h浙江大学研究生学位课程546.2.2 Adams四步四阶四步四阶预估校正算法预估校正算法浙江大学研究生学位课程55 6.2.3 Runge-Kutta 方法方法浙江大学研究生学位课程566.2.3 浙江大学研究生学位课程576.2.3 浙江大学研究生学位课程586.2.3 浙江大学研究生学位课程596.2.3 六个未知数，二个自由六个未知数，二个自由，故可取故可取故：故：浙江大学研究生学位课程606.2.3 N=4:四级四阶四级四阶R-K方法方法 -最常用的古典最常用的古典Runge-Kutta方法。方法。增加计算函数值的次数（级）与提高精增加计算函数值的次数（级）与提高精 度（阶）的关系见下表：度（阶）的关系见下表：表表 6.7 Runge-Kutta方法中方法中级级与与阶阶的关系的关系浙江大学研究生学位课程616.2.3 表表 6.8 各种解法在例题中的结果比较各种解法在例题中的结果比较浙江大学研究生学位课程626.2.3 （2）.单步法，自动起步单步法，自动起步 （3）.易改为变步长易改为变步长 （4）.绝对稳定区域较同阶线性多步法大绝对稳定区域较同阶线性多步法大 （5）.计算工作量较大，有时大于隐式方法计算工作量较大，有时大于隐式方法 （6）.估计误差不易估计误差不易 绝对稳定性讨论绝对稳定性讨论 浙江大学研究生学位课程636.2.3 表表 6.9 各级各级R-K方法的绝对稳定区域方法的绝对稳定区域浙江大学研究生学位课程646.2.3 表表 6.10 不同步长对精度的影响不同步长对精度的影响浙江大学研究生学位课程656.2.3 浙江大学研究生学位课程666.2.3 浙江大学研究生学位课程67 P阶图图 6.10 变步长Runge-Kutta方法框图6.2.3 浙江大学研究生学位课程686.2.3 浙江大学研究生学位课程696.2.4 出发值的计算出发值的计算浙江大学研究生学位课程706.2.4 浙江大学研究生学位课程71质量控制质量控制 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.5 HCALL 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误差计算一个单步计算两个半步长计算置初始值111126.2.4 浙江大学研究生学位课程72PARAMETER(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)RETURNEND6.2.4 浙江大学研究生学位课程73IF(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 步进程序（续）步进程序（续）136.2.4 浙江大学研究生学位课程74 计算结果计算结果 NOK=30 NBAD=0 KOUNT=156.2.4 浙江大学研究生学位课程75四阶四阶 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)CONTINUERETURNEND141312116.2.4 浙江大学研究生学位课程76用四阶用四阶 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结束返回结束返回12136.2.4 浙江大学研究生学位课程77带自适应步长控制的四阶带自适应步长控制的四阶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存贮中间结果存贮中间结果131211置初值6.2.4 浙江大学研究生学位课程78IF (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 步数太多，越界步数太多，越界RETURNEND1615146.2.4 浙江大学研究生学位课程796.2.4 浙江大学研究生学位课程80 6.2.5 方程组与刚性问题方程组与刚性问题 （Stiff Sets of Equation）1.考虑考虑2阶常微分方程组的初值问题：阶常微分方程组的初值问题：解此方程组的解此方程组的Euler方法方法浙江大学研究生学位课程816.2.5 浙江大学研究生学位课程826.2.5 可能是复数可能是复数浙江大学研究生学位课程836.2.5 绝对收敛绝对收敛。Adams内插法的绝对稳定区域内插法的绝对稳定区域：图图 6.11 Adams内插法的绝对稳定区域内插法的绝对稳定区域浙江大学研究生学位课程846.2.5图图6.12 Adams外推法的绝对稳定区域外推法的绝对稳定区域图图 6.13 R-K法的绝对稳定区域法的绝对稳定区域浙江大学研究生学位课程856.2.5 2.刚性方程组刚性方程组（Stiff Equations）浙江大学研究生学位课程866.2.5 浙江大学研究生学位课程876.2.5 定义定义：刚性方程刚性方程：一般坏条件问题一般坏条件问题严重坏条件问题严重坏条件问题刚性比刚性比：浙江大学研究生学位课程886.2.5 图图 6.14 A-稳定区域稳定区域浙江大学研究生学位课程896.2.5 图图 6.15浙江大学研究生学位课程906.2.5浙江大学研究生学位课程916.2.5 浙江大学研究生学位课程926.2.5 浙江大学研究生学位课程936.2.4 浙江大学研究生学位课程94Euler法一阶中点法二阶(Runge-Kutta)四级四阶RungeKuttaEuler 预估-校正法二阶6.2.4 浙江大学研究生学位课程956.3 边值问题解法边值问题解法浙江大学研究生学位课程96 6.3.1 试射法试射法（Shooting）图图 6.16 试射法求解示意试射法求解示意浙江大学研究生学位课程976.3.1浙江大学研究生学位课程986.3.1浙江大学研究生学位课程99 6.3.2 差分方法差分方法（Difference Methods）(Relaxation Methods)True solution2nd iteration1st iterationInitial guess图图 6.17 差分方法示意差分方法示意浙江大学研究生学位课程1006.3.2 浙江大学研究生学位课程1016.3.2浙江大学研究生学位课程1026.3.2图图 6.18 修正值矩阵修正值矩阵浙江大学研究生学位课程103 块对角线性代数方程组的快速求解块对角线性代数方程组的快速求解 利用矩阵的特殊结构使总运算次数达到极小，利用矩阵的特殊结构使总运算次数达到极小，并使矩阵系数的存贮达到极小并使矩阵系数的存贮达到极小图图 6.19 Gauss 消去法的目标结构消去法的目标结构6.3.2浙江大学研究生学位课程104D-有待于对角化有待于对角化A-将要改动将要改动S-需要存贮需要存贮Z-将要变为将要变为0 a)b)a)b)6.3.2图图 6.20 特殊矩阵结构特殊矩阵结构 1 浙江大学研究生学位课程105 a)b)包括四步运算包括四步运算1.使用前一步的结果，简化部分元素为零使用前一步的结果，简化部分元素为零2.消去剩下子块正方结构为单位阵消去剩下子块正方结构为单位阵3.存储以后要使用的非零系数存储以后要使用的非零系数4.回代回代6.3.2图图 6.22 特殊矩阵结构特殊矩阵结构 2浙江大学研究生学位课程106 6.3.3 样条函数在两点边值问题中的应用样条函数在两点边值问题中的应用 考察考察 二阶线性微分方程二阶线性微分方程 的边值问题的边值问题浙江大学研究生学位课程1076.3.3浙江大学研究生学位课程1086.3.3 浙江大学研究生学位课程1096.3.3 浙江大学研究生学位课程1106.3.3浙江大学研究生学位课程1116.3.3浙江大学研究生学位课程112 相平面上的轨迹为相平面上的轨迹为1914年-1923年捕获的鱼中，食用鱼所占的比例(%)图图 6.23 相平面上的轨迹图相平面上的轨迹图6.3.3浙江大学研究生学位课程1136.3.3浙江大学研究生学位课程114公式简单，但精度不好公式简单，但精度不好编程简单，使用方便编程简单，使用方便计算量少，但稳定性差计算量少，但稳定性差计算量大，但稳定性好计算量大，但稳定性好6.3.3总总 结结浙江大学研究生学位课程115第六章第六章 习题习题浙江大学研究生学位课程116