浙江大学研究生学位课程,实用数值计算方法,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, YTEM