《常微分方程的数值解法》由会员分享，可在线阅读，更多相关《常微分方程的数值解法（93页珍藏版）》请在金锄头文库上搜索。
1、Tel： 86613747E-mail： 授课: 68学分：4第七章 常微分方程的数值解法 问题提出倒葫芦形状容器壁上的刻度问题.对于如图所示圆 柱形状容器壁上的容积刻度,可以利用圆柱体体积公式其中直径D为常数.由于体积V与相对于容器底部的 任意高度H的函数关系明确,因此在容器上可以方便 地标出容器刻度,而对于几何形状不是规则的容器, 比如倒葫芦形状容器壁上如何标出刻度呢?下表是 经过测量得到部分容器高度与直径的关系.x1o根据上表的数据,可以拟合出倒 葫芦形状容器的图,建立如图所 示的坐标轴后,问题即为如何根 据任意高度x标出容器体积V的 刻度,由微元思想分析可知H 0 0.2 0.4 0.
2、6 0.8 1.0D 0.04 0.11 0.26 0.56 1.04 1.17其中x表示高度,直径D是高度x的函数,记为D(x), 因此得到如下微分方程初值问题第七章 常微分方程的数值解法 只要求解上述方程,就可求出体积V与高度x之间的 函数关系,从而可标出容器壁上容积的刻度,但问题 是(7.0)中的函数D(x)无解析表达式,我们根本无法 求出其解析解.怎样求解不能用解析方法求解的微分方程初制 值问题,就是本章要讨论的重点.( 7.0 )第七章 常微分方程的数值解法 7.1 引言包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中, 自变量的个数只有一个, 称为常微分方
3、程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数都是一次的,则称它是线性的,否则称为非线性的。 在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。 譬如 这个一阶微分方程就不能用初等函数及其积分来表达它的解。 再如，方程 的解 ,虽然有表可查,但对于表上没有给出 的值,仍需插值方法来计算从实际问题当中归纳出来的微分方程，通常主要依靠数值
4、解法来解决。本章主要讨论一阶常微分方程初值问题 ( 7.1 ) 在区间a x b上的数值解法。 可以证明,如果函数在带形区域 R=axb,-y内连续，且关于y满足李普希兹 (Lipschitz)条件，即存在常数L(它与x,y无关)使 对R内任意两个 都成立,则方程( 7.1 )的解在a, b上存在且唯一。 常微分方程表示方法在微分方程中, 自变量的个数只有一个, 称为常微分方程.自变量的个数为两个或两个以上的微分方程叫偏微分方程。例如当x=0时,y=1,可得c=1特解当x=0时,y=1,可得c=-1特解两边积分通解常微分方程解法回顾-例：设降落伞从跳伞塔下落后，所受空气阻力与速度 成正比，并设
5、降落伞离开跳伞塔时（t=0）速度为零，求降落伞下落速度与时间的函数关系。 解：设降落伞下落速度为v(t),降落伞在空中下落时， 同时受到重力P与阻力R的作用，重力大小为mg，方向 与v一致；阻力大小为kv（k为比例系数），方向与v相 反，从而降落伞所受外力为F=mg-kv 根据牛顿第二定律 F=ma (a 为加速度），得函数v(t)应满足的方程为按题意，初始条件为分离变量后得两端积分考虑到mg-kv0即或将初始条件代如入上式得于是所求特解为7.2 数值方法的基本思想对常微分方程初值问题(7.1)式的数值解法，就是要算出精确解y(x)在区间a,b上的一系列离散节点 处的函数值 的近似值。相邻两个
6、节点的间距 称为步长，步长可以相等，也可以不等。本章总是假定h为定数，称为定步长，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。 对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点，它们都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息 计算 的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题中的导数 进行不同的离散化处理。 对于初值问题的数值解法，首先要解决的问题就是如何对微分方程进行离散化，建立求数值解的递推公式。递推公式通 常有两类，
7、一类是计算yi+1时只用到xi+1, xi 和yi， 即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法；其代表是龙格库塔法。另一类是计算yi+1时，除用到xi+1,xi和yi以外，还要用到 ，即前面k步的值，此类方法称为多步法；其代表是亚当斯法。 7.3 欧拉（Euler）法 7.3.1 Euler公式欧拉（Euler）方法是解初值问题的最简单的数值方法。初值问题的解y=y(x)代表通过点 的一条称之为 微分方程的积分曲线。积分曲线上每一点的切线的斜率 等于函数 在 这点的值。 Euler法的求解过程是:从初始点P0(即点(x0,y0)出发,作积分曲线y=y(x)在P0点上
8、切线 (其斜率为),与x=x1直线相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值,如上图所示。过点(x0,y0),以f(x0,y0)为斜率的切线 方程为 当 时,得 这样就获得了P1点的坐标。 同样, 过点P1(x1,y1),作积分曲线y=y(x)的切线交直线x=x2于P2点,切线 的斜率 =直线方程为当 时,得 当 时,得由此获得了P2的坐标。重复以上过程,就可获得一系列的点:P1,P1,Pn。对已求得点以 = 为斜率作直线 取从图形上看,就获得了一条近似于曲线y=y(x)的折线 。这样,从x0逐个算出对应的数值解 通常取 (常数),则Euler法的计算格式 i=0,1,n
9、 ( 7.2 ) 还可用数值微分、数值积分法和泰勒展开法推导Euler格式。以数值积分为例进行推导。将方程 的两端在区间 上积分得， 选择不同的计算方法计算上式的积分项,就会得到不同的计算公式。 (7.3 )公式推导用左矩形方法计算积分项 代入(7.3)式,并用yi近似代替式中y(xi)即可得到向前欧拉（Euler）公式 由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。 例7.1 用欧拉法解初值问题 取步长h=0.2 ,计算过程保留4位小数 解: h=0.2, 欧拉迭代格式 当 k=0, x1=0.2时，已知x0=0，y0=1，有y(0.2)y1=0.21(401)0.8
10、 当 k=1, x2=0.4时，已知x1 =0.2， y1 =0.8，有y(0.4) y2 =0.20.8(40.20.8)0.6144 当 k=2, x3 =0.6时，已知x2 =0.4， y2 =0.6144，有y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 例7.2 用Euler法求初值问题y= x y 2y (0) = 0 的数值解 （ 取 h = 0.1 n = 5 )解： f (x,y) = x - y 2 ; x0 = y0= 0 ; h=0.1 由Euler法的递推公式得：yn+1 = yn + 0.1 (xn y2n )yn = 0 n = 0
11、, 1, 2, 3, 4, 5由上式计算所得数据n 0 1 2 3 4 5 xn 0 0.1 0.2 0.3 0.4 0.5 yn 0 0 0.01 0.02999 0.05999 0.0995 二、隐式Euler ( 尤拉)格式7.3.2 梯形公式为了提高精度,对方程 的两端在区间上积分得，改用梯形方法计算其积分项，即 ( 7.4 ) 代入(7.4)式,并用近似代替式中即可得到梯形公式 ( 7.5 ) 由于数值积分的梯形公式比矩形公式的精度高， 因此梯形公式（7.5）比欧拉公式( 7.2 )的精度高一 个数值方法。 ( 7.5 ) (7.5)式的右端含有未知的yi+1,它是一个关于yi+1的
12、函数方程,这类数值方法称为隐式方法。相反地,欧拉法是关于yi+1的一个直接的计算公式， 这类数值方法称为显式方法。 例7.3 用梯形公式求下面初值问题的解在x=0.01上的值y(0.01)y= y y (0) = 1 解：取h=0.01, 由梯形公式得 y=ex e0.01=1.0100501677.3.3 两步欧拉公式对方程 的两端在区间上 积分得 ( 7.6 ) 改用中矩形公式计算其积分项，即 代入上式,并用yi近似代替式中y(xi)即可得到两步欧拉公式 ( 7.7 ) 前面介绍过的数值方法,无论是欧拉方法,还是梯形方法，它们都是单步法,其特点是在计算yi+1时只用到前一步的信息yi;可是
13、公式(7.7)中除了yi外,还用到更前一步的信息yi-1,即调用了前两步的信息,故称其为两步欧拉公式 7.3.4 欧拉法的局部截断误差衡量求解公式好坏的一个主要标准是求解公式的精度, 因此引入局部截断误差和阶数的概念。定义7.1 在yi准确的前提下, 即 时, 用数值方法计算yi+1的误差 , 称为该数值方法计算时yi+1的局部截断误差。对于欧拉公式，假定 ，则有而将真解y(x)在xi处按二阶泰勒展开 因此有 定义7.2 数值方法的局部截断误差为 ,则称这种数值方法的阶数是P。步长(h N 结束。 （2）改进欧拉法的流程图 (3) 程序实现(见附录A A-15 改进欧拉法计算常微分方程初值问题
14、 ) 例7.2 用改进欧拉法解初值问题 区间为0,1,取步长h=0.1 解: 改进欧拉法的具体形式 本题的精确解为 ,计算见P158列表所示 例7.3 对初值问题 证明用梯形公式求得的近似解为 并证明当步长h0时,yn收敛于精确解 证明: 解初值问题的梯形公式为 整理成显式 反复迭代,得到 由于 ，有 证毕 7.4 龙格-库塔（Runge-Kutta）法 7.4.1 龙格-库塔(Runge-Kutta)法的基本思想Euler公式可改写成 则yi+1的表达式与 y(xi+1)的Taylor展开式的前两项 完全相同,即局部截断误差为 。改进的Euler公式又可改写成 上述两组公式在形式上有一个共同点:都是用f(x,y)在某些点上值的线性组合得出y(xi+1)的近似值yi+1,而且增加计算的次数f(x,y)的次数,可提高截断误
