首页 数值分析-第七章常微分方程数值解

数值分析-第七章常微分方程数值解

举报
开通vip

数值分析-第七章常微分方程数值解nullnullTel: 86613747 E-mail: lss@zjtcm.net 授课: 68 学分:4null第七章 常微分方程的数值解法 问题提出 倒葫芦形状容器壁上的刻度问题.对于如图所示圆柱形状容器壁上的容积刻度,可以利用圆柱体体积公式其中直径D为常数.由于体积V与相对于容器底部的任意高度H的函数关系明确,因此在容器上可以方便地标出容器刻度,而对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢?下表是经过测量得到部分容器高度与直径的关系.nullx1o根据上表的数据,可以...

数值分析-第七章常微分方程数值解
nullnullTel: 86613747 E-mail: lss@zjtcm.net 授课: 68 学分:4null第七章 常微分方程的数值解法 问题提出 倒葫芦形状容器壁上的刻度问题.对于如图所示圆柱形状容器壁上的容积刻度,可以利用圆柱体体积公式其中直径D为常数.由于体积V与相对于容器底部的任意高度H的函数关系明确,因此在容器上可以方便地标出容器刻度,而对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢?下 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 是经过测量得到部分容器高度与直径的关系.nullx1o根据上表的数据,可以拟合出倒葫芦形状容器的图,建立如图所示的坐标轴后,问题即为如何根据任意高度x标出容器体积V的刻度,由微元思想分析可知H 0 0.2 0.4 0.6 0.8 1.0 D 0.04 0.11 0.26 0.56 1.04 1.17其中x表示高度,直径D是高度x的函数,记为D(x),因此得到如下微分方程初值问题null第七章 常微分方程的数值解法 只要求解上述方程,就可求出体积V与高度x之间的函数关系,从而可标出容器壁上容积的刻度,但问题是(7.0)中的函数D(x)无解析表达式,我们根本无法求出其解析解. 怎样求解不能用解析方法求解的微分方程初制值问题,就是本章要讨论的重点.( 7.0 )null第七章 常微分方程的数值解法 7.1 引言 包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中, 自变量的个数只有一个, 称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数 都是一次的,则称它是线性的,否则称为非线性的。 null 在高等 数学 数学高考答题卡模板高考数学答题卡模板三年级数学混合运算测试卷数学作业设计案例新人教版八年级上数学教学计划 中,对于常微分方程的求解,给出了一些典型方程求解析解的基本方法,如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是有限的,大多数的常微分方程是不可能给出解析解。 譬如 这个一阶微分方程就不能用初等函数及其积分来表达它的解。 null再如,方程 的解 ,虽然有表可查,但对于表上没有给出 的值,仍需插值方法来计算 null从实际问题当中归纳出来的微分方程,通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题 ( 7.1 ) 在区间a ≤ x ≤ b上的数值解法。 可以 证明 住所证明下载场所使用证明下载诊断证明下载住所证明下载爱问住所证明下载爱问 ,如果函数在带形区域 R=a≤x≤b, -∞<y<∞}内连续,且关于y满足李普希兹(Lipschitz)条件,即存在常数L(它与x,y无关)使 对R内任意两个 都成立,则方程( 7.1 )的解 在a, b上存在且唯一。 null常微分方程表示方法null在微分方程中, 自变量的个数只有一个, 称为 常微分方程. 自变量的个数为两个或两个以上的微分方程 叫偏微分方程。null例如null常微分方程解法回顾当x=0时,y=1,可得c=1特解当x=0时,y=1,可得c=-1特解两边积分通解-null例:设降落伞从跳伞塔下落后,所受空气阻力与速度 成正比,并设降落伞离开跳伞塔时(t=0)速度为 零,求降落伞下落速度与时间的函数关系。 解:设降落伞下落速度为v(t),降落伞在空中下落时,同时受到重力P与阻力R的作用,重力大小为mg,方向与v一致;阻力大小为kv(k为比例系数),方向与v相反,从而降落伞所受外力为 F=mg-kv 根据牛顿第二定律 F=ma (a 为加速度),得函数v(t)应满足的方程为 按题意,初始条件为分离变量后得两端积分考虑到mg-kv>0即或将初始条件代如入上式得于是所求特解为nullnull7.2 数值方法的基本思想 对常微分方程初值问题(7.1)式的数值解法,就是要算出精确解y(x)在区间a,b上的一系列离散节点 处的函数值 的近似值 。相邻两个节点的间距 称为步长,步 长可以相等,也可以不等。本章总是假定h为定数,称为定步长,这时节点可表示为 数值解法需要把连续性的问题加以离散化,从而求出离散节点的数值解。 null 对常微分方程数值解法的基本出发点就是离散化。其数值解法有两个基本特点,它们都采用“步进式”,即求解过程顺着节点排列的次序一步一步地向前推进,描述这类算法,要求给出用已知信息 计算 的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法,对初值问题 中的导数 进行不同的离散化处理。 null对于初值问题 的数值解法,首先要解决的问题就是如何对微分方程进行离散化,建立求数值解的递推公式。递推公式通常有两类,一类是计算yi+1时只用到xi+1, xi 和yi, 即前一步的值,因此有了初值以后就可以逐步往下计算,此类方法称为单步法;其代表是龙格—库塔法。另一类是计算yi+1时,除用到xi+1,xi和yi以外,还要用到 ,即前面k步的值,此类方法称为多步法;其代表 是亚当斯法。 null7.3 欧拉(Euler)法 7.3.1 Euler公式 欧拉(Euler)方法是解初值问题的最简单的数值方法。初值问题 的解y=y(x)代表通过点 的一条称之为微分方程的积分曲线。积分曲线上每一点 的切线的斜率 等于函数 在这点的值。 nullEuler法的求解过程是:从初始点P0(即点(x0,y0))出发, 作积分曲线y=y(x)在P0点上切线 (其斜率为 ),与x=x1直线相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值,如上图所示。过点(x0,y0),以f(x0,y0)为斜率的切线方程为 当 时,得 这样就获得了P1点的坐标。 null同样, 过点P1(x1,y1),作积分曲线y=y(x)的切线 交直线x=x2于P2点,切线 的斜率 = 直线方程为当 时,得 null当 时,得 由此获得了P2的坐标。重复以上过程,就可获得一系列的点:P1,P1,…,Pn。对已求得点 以 = 为斜率作直线 取null 从图形上看,就获得了一条近似于曲线y=y(x) 的折线 。这样,从x0逐个算出 对应的数值解 null通常取 (常数),则Euler法的计算格式 i=0,1,…,n ( 7.2 ) 还可用数值微分、数值积分法和泰勒展开法推导Euler格式。以数值积分为例进行推导。 将方程 的两端在区间 上积分得, 选择不同的计算方法计算上式的积分项 ,就会得到不同的计算公式。 (7.3)null 公式推导null null 用左矩形方法计算积分项 代入(7.3)式,并用yi近似代替式中y(xi)即可得到向前欧拉(Euler)公式 由于数值积分的矩形方法精度很低,所以欧拉(Euler)公式当然很粗糙。 null null例7.1 用欧拉法解初值问题 取步长h=0.2 ,计算过程保留4位小数 解: h=0.2, 欧拉迭代格式 当 k=0, x1=0.2时,已知x0=0,y0=1,有 y(0.2)y1=0.2×1(4-0×1)=0.8 当 k=1, x2=0.4时,已知x1 =0.2, y1 =0.8,有 y(0.4) y2 =0.2×0.8×(4-0.2×0.8)=0.6144 当 k=2, x3 =0.6时,已知x2 =0.4, y2 =0.6144,有 y(0.6) y3=0.2×0.6144×(4-0.4×0.6144)=0.4613 null例7.2 用Euler法求初值问题 y'= x – y 2 y (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, 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 … null二、隐式Euler ( 尤拉)格式null7.3.2 梯形公式 为了提高精度,对方程 的两端在区间上 积分得, 改用梯形方法计算其积分项,即 ( 7.4 ) 代入(7.4)式,并用近似代替式中即可得到梯形公式 ( 7.5 ) 由于数值积分的梯形公式比矩形公式的精度高,因此梯形公式(7.5)比欧拉公式( 7.2 )的精度高一个数值方法。 null( 7.5 ) (7.5)式的右端含有未知的yi+1,它是一个关于yi+1的函数方程,这类数值方法称为隐式方法。相反地,欧拉法是关于yi+1的一个直接的计算公式, 这类数值方法称为显式方法。 null例7.3 用梯形公式求下面初值问题的解 在x=0.01上的值y(0.01) y'= y y (0) = 1 解:取h=0.01, 由梯形公式得 y=ex e0.01=1.010050167null7.3.3 两步欧拉公式 对方程 的两端在区间上 积分得 ( 7.6 ) 改用中矩形公式计算其积分项,即 代入上式,并用yi近似代替式中y(xi)即可得到两步欧拉公式 ( 7.7 ) null 前面介绍过的数值方法,无论是欧拉方法,还是梯形方法,它们都是单步法,其特点是在计算yi+1时只用到前一步的信息yi;可是公式(7.7)中除了yi外,还用到更前一步的信息yi-1,即调用了前两步的信息,故称其为两步欧拉公式 null7.3.4. 欧拉法的局部截断误差 衡量求解公式好坏的一个主要 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 是求解公式的精度, 因此引入局部截断误差和阶数的概念。 定义7.1 在yi准确的前提下, 即 时, 用数值方法计算yi+1的误差 , 称为该数值方法计算时yi+1的局部截断误差。 对于欧拉公式,假定 ,则有 而将真解y(x)在xi处按二阶泰勒展开 因此有 null定义7.2 数值方法的局部截断误差为 ,则称这种数值方法的阶数是P。步长(h<1) 越小,P越高, 则局部截断误差越小,计算精度越高。欧拉公式的局部截断误差为 , 欧拉方法仅为一阶方法。 两步欧拉公式比欧拉公式精度也是高一个数值方法,设 , 前两步准确,则两步欧拉公式 把y(xi-1)在xi处展开成Taylor级数,即 将在xn处按二阶泰勒展开略去余项用近似值yn代替y(xn)有null由 再将y(xi+1)在xi处进行泰勒展开(7.8)(7.9)所以, 由(7.8)和(7.9)可得两步欧拉公式的局部截断 误差为, 即 null7.3.5 改进的欧拉公式 显式欧拉公式计算工作量小,但精度低。梯形公式虽提高了精度,但为隐式公式,需用迭代法求解,计算工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。 先用欧拉公式(7.2)求出一个初步的近似值 ,称为预测值, 它的精度不高, 再用梯形公式(7.5)对它校正一次,即迭代一次,求得yi+1,称为校正值, 这种预测-校正方法称为改进的欧拉公式:(7.10) 预测 校正 null 可以证明,公式(7.10)的精度为二阶。这是一种 一步显式格式,它可以表示为嵌套形式。 (7.11)或者表示成下列平均化形式 (7.12) null7.3.6 改进欧拉法算法实现 (1)计算步骤 ① 输入 , h , N ② 使用以下改进欧拉法公式进行计算 ③ 输出 ,并使 转到 ② 直至n > N 结束。 null(2)改进欧拉法的流程图 null(3) 程序实现(见附录A A-15 改进欧拉法计算常微 分方程初值问题 ) 例7.2 用改进欧拉法解初值问题 区间为0,1,取步长h=0.1 解: 改进欧拉法的具体形式 本题的精确解为 ,计算见P158列表所示 null例7.3 对初值问题 证明用梯形公式求得的近似解为 并证明当步长h0时,yn收敛于精确解 证明: 解初值问题的梯形公式为 ∵ ∴ 整理成显式 反复迭代,得到 ∵ ∴ null由于 ,有 ∴ 证毕 null7.4 龙格-库塔(Runge-Kutta)法 7.4.1 龙格-库塔(Runge-Kutta)法的基本思想 Euler公式可改写成 则yi+1的表达式与 y(xi+1)的Taylor展开式的前两项完全相同,即局部截断误差为 。 改进的Euler公式又可改写成 null 上述两组公式在形式上有一个共同点:都是用f(x,y)在某些点上值的线性组合得出y(xi+1)的近似值yi+1,而且增加计算的次数f(x,y)的次数,可提高截断误差的阶。如欧拉公式:每步计算一次f(x,y)的值,为一阶方法。改进欧拉公式需计算两次f(x,y)的值,它是二阶方法。它的局部截断误差为 。 null 于是可考虑用函数f(x,y)在若干点上的函数值的线性组合来构造近似公式,构造时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合,从而使近似公式达到所需要的阶数。既避免求偏导,又提高了计算方法精度的阶数。或者说,在 这一步内多预报几个点的斜率值,然后将其加权平均作为平均斜率,则可构造出更高精度的计算格式,这就是龙格—库塔(Runge-Kutta)法的基本思想。 null7.4.2 二阶龙格—库塔法 在 上取两点xi和 ,以该两点处的斜率值k1和k2的加权平均(或称为线性组合)来求取平均斜率k*的近似值K,即 式中:k1为xi点处的切线斜率值, k2为 点处的切线斜率值,比照改进的欧拉法,将 视为 ,即可得 对常微分方程初值问题(7.1)式的解 y=y(x),根据微分中值定理,存在点 ,使得 null式中 K可看作是y=y(x)在区间 上的平均斜率。所以可得计算公式为: (7.14) 将y(xi)在x=xi处进行二阶Taylor展开: (7.15) 也即 (7.13)null将 在x=xi处进行一阶Taylor展开: 将以上结果代入(7.14)得: (7.16) 对式(7.15)和(7.16)进行比较系数后可知,只要 (7.17) 成立,格式(7.14)的局部截断误差就等于有2阶 精度null式(7.17)中具有三个未知量,但只有两个方程,因而有无穷多解。若取 ,则p=1,这是无穷多解中的一个解,将以上所解的值代入式(7.14)并改写可得 不难发现,上面的格式就是改进的欧拉格式。凡满足条件式(7.17)有一簇形如上式的计算格式,这些格式统称为二阶龙格—库塔格式。因此改进的欧拉格式是众多的二阶龙格—库塔法中的一种特殊格式。 null若取 ,则 ,此时二阶龙格-库塔 法的计算公式为 此计算公式称为变形的二阶龙格—库塔法。式中 为区间 的中点。 null7.4.3 三阶龙格-库塔法 为了进一步提高精度,设除 外再增加一点 并用三个点 , , 的斜率k1,k2,k3加权平均 得出平均斜率k*的近似值,这时计算格式具有形式: (7.18) 为了预报点 的斜率值k3,在区间 内有两 个斜率值k1和k2可以用,可将k1,k2加权平均得出 上的平均斜率,从而得到 的预报值 null于是可得 运用Taylor展开方法选择参数 ,可以使格式(7.18)的局部截断误差为 ,即具有三阶精度,这类格式统称为三阶龙格—库塔方法。下列是其中的一种,称为库塔(Kutta)公式。 (7.19) null7.4.4 四阶龙格—库塔法 如果需要再提高精度,用类似上述的处理方法,只需在区间 上用四个点处的斜率加权平均作为平均斜率k*的近似值,构成一系列四阶龙格—库塔公式。具有四阶精度,即局部截断误差是 。 由于推导复杂,这里从略,只介绍最常用的一种四阶经典龙格—库塔公式。 (7.20) null7.4.5 四阶龙格—库塔法算法实现 (1) 计算步骤 ① 输入 ,h,N ② 使用龙格—库塔公式(7.20)计算出y1 ③ 输出 ,并使 转到 ② 直至n > N 结束。 null(2) 四阶龙格—库塔算法流程图null 程序实现(见附录A A-16 四阶龙格-库塔法计 算常微分方程初值问题) 例7.4 取步长h=0.2,用经典格式求解初值问题 解: 由四阶龙格-库塔公式可得 null可同样进行其余yi的计算。本例方程的解为 ,数值解yi与准确解y(xi)的对照表见教材P163所示,从表中看到所求的数值解具有4位有效数字。 龙格—库塔方法的推导基于Taylor展开方法,因而它要求所求的解具有较好的光滑性。如果解的光滑性差,那么,使用四阶龙格—库塔方法求得的数值解,其精度可能反而不如改进的欧拉方法。在实际计算时,应当针对问题的具体特点选择合适的算法。 null7.4.6 变步长的龙格-库塔法 在微分方程的数值解中,选择适当的步长是非常重要的。单从每一步看,步长越小,截断误差就越小;但随着步长的缩小,在一定的求解区间内所要完成的步数就增加了。这样会引起计算量的增大,并且会引起舍入误差的大量积累与传播。因此微分方程数值解法也有选择步长的问题。 以经典的四阶龙格-库塔法(7.20)为例。从节点xi出发,先以h为步长求出一个近似值,记为 ,由于局部截断误差为 ,故有 当h值不大时,式中的系数c可近似地看作为常数。null然后将步长折半,即以为 步长,从节点xi出发,跨两步到节点xi+1,再求得一个近似值 ,每跨一步的截断误差是 ,因此有这样 由此可得 这表明以 作为 的近似值,其误差可用步长折半前后两次计算结果的偏差 来判断所选步长是否适当null当要求的数值精度为ε时: (1)如果Δ>ε,反复将步长折半进行计算,直至Δ<ε为止,并取其最后一次步长的计算结果作为 (2)如果Δ<ε,反复将步长加倍,直到Δ>ε为止,并以上一次步长的计算结果作为 。 这种通过步长加倍或折半来处理步长的方法称为变步长法。表面上看,为了选择步长,每一步都要反复判断Δ,增加了计算工作量,但在方程的解y(x)变化剧烈的情况下,总的计算工作量得到减少,结果还是合算的。 null7.5 亚当姆斯方法 7.5.1 亚当姆斯格式 龙格-库塔方法是一类重要算法,但这类算法在每一步都需要先预报几个点上的斜率值,计算量比较大。考虑到计算yi+1之前已得出一系列节点上 的斜率值,能否利用这些已知值来减少计算量呢? 这就是亚当姆斯(Adams)方法的 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 思想。 null 设用xi,xi+1两点的斜率值加权平均作为区间 上的平均斜率,有计算格式 (7.21) 选取参数λ,使格式(7.21)具有二阶精度。 null将 在xi处Taylor展开 代入计算格式(7.21)化简,并假设, 因此有 与y(xi+1)在xi处的Taylor展开式 相比较, 需取 才使格式(7.21)具有二阶精度。这样导出的计算格式 称之为二阶亚当姆斯格式。类似地可以导出三阶亚当姆斯格式。 null和四阶亚当姆斯格式。 (7.22) 这里和下面均记 上述几种亚当姆斯格式都是显式的,算法比较简单,但用节点 的斜率值来预报区间 上的平均斜率是个外推过程,效果不够理想。为了进一步改善精度,变外推为内插,即增加节点xi+1的斜率值来得出 上的平均斜率。譬如考察形如 null(7.23) 的隐式格式,设(7.23)右端的 Taylor展开有 可见要使格式(7.23)具有二阶精度,需令 , 这样就可构造二阶隐式亚当姆斯格式 其实是梯形格式。类似可导出三阶隐式亚当姆斯格式 null和四阶隐式亚当姆斯格式 (7.24) 7.5.2 亚当姆斯预报-校正格式 参照改进的欧拉格式的构造方法,以四阶亚当姆斯为例,将显式(7.22)和隐式(7.24)相结合,用显式公式做预报,再用隐式公式做校正,可构成亚当姆斯预报-校正格式 (7.25) 预报: 校正: null 这种预报-校正格式是四步法,它在计算yi+1时不但用到前一步的信息 ,而且要用到再前面三步的信息 ,因此它不能自行启动。在实际计算时,可借助于某种单步法,譬如四阶龙格—库塔法提供开始值 。 null例7.5 取步长h=0.1,用亚当姆斯预报-校正公式求解 初值问题 的数值解。 解: 用四阶龙格-库塔公式求出发值 ,计算得: 表中的 ,yi和y(xi)分别为预报值、校正值和准确解( ),以比较计算结果的精度。 再使用亚当姆斯预报-校正公式(7.25),见教材P166列表算得其余的计算结果null 7.6 算法的稳定性及收敛性 7.6.1 稳定性 稳定性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是用差分格式进行计算的,而在差分方程的求解过程中,存在着各种计算误差,这些计算误差如舍入误差等引起的扰动,在传播过程中,可能会大量积累,对计算结果的准确性将产生影响。这就涉及到算法稳定性问题。 null 当在某节点上xi的yi值有大小为δ的扰动时,如果在其后的各节点 上的值yi产生的偏差都不大于δ,则称这种方法是稳定的。 稳定性不仅与算法有关,而且与方程中函数f(x,y)也有关,讨论起来比较复杂。为简单起见,通常只针对模型方程来讨论。一般方程若局部线性化,也可化为上述形式。模型方程相对比较简单,若一个数值方法对模型方程是稳定的,并不能保证该方法对任何方程都稳定,但若某方法对模型方程都不稳定,也就很难用于其他方程的求解。null先考察显式Euler方法的稳定性。模型方程 的Euler公式为 将上式反复递推后,可得 或式中null 要使yi有界,其充要条件为 即 由于 ,故有 (7.26) 可见,如欲保证算法的稳定,显式Euler格式的步长h的选取要受到式(7.26)的限制。 的绝对值越大,则限制的h值就越小。 用隐式Euler格式,对模型方程 的计算公式为,可化为 null由于 ,则恒有 ,故恒有 因此,隐式Euler格式是绝对稳定的(无条件稳定)的(对任何h>0)。 7.6.2 收敛性 常微分方程初值问题的求解,是将微分方程转化为差分方程来求解,并用计算值yi来近似替代y(xi),这种近似替代是否合理,还须看分割区间 的长度h越来越小时,即 时, 是否成立。若成立,则称该方法是收敛的,否则称为不收敛。 这里仍以Euler方法为例,来分析其收敛性。 Euler格式null设 表示取 时, 按Euler公式的计算结果, 即 Euler方法局部截断误差为:设有常数 ,则 (7.27) 总体截断误差 (7.28) 又 由于f(x,y)关于y满足李普希茨条件,即 null代入(7.28)上式,有 再利用式(7.27),式(7.28) 即 上式反复递推后,可得 (7.29) null设 (T为常数) 因为 所以 把上式代入式(7.29),得 若不计初值误差,即 ,则有 (7.30) 式(7.30)说明,当时 , ,即 ,所以Euler方法是收敛的,且其收敛速度为 ,即具有一阶收敛速度。同时还说明Euler方法的整体截断误差为 ,因此算法的精度为一阶。 null7.7 一阶常微分方程组与高阶方程 我们已介绍了一阶常微分方程初值问题的各种数值解法,对于一阶常微分方程组,可类似得到各种解法,而高阶常微分方程可转化为一阶常微分方程组来求解。  7.7.1 一阶常微分方程组 对于一阶常微分方程组的初值问题 (7.31) 可以把单个方程 中的f 和y看作向量来处理,这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。 null设 为节点上的近似解, 则有改进的Euler格式为 预报:校正: (7.32) 又,相应的四阶龙格—库塔格式(经典格式)为 (7.33) null式中 (7.34) null 把节点xi上的yi和zi值代入式(7.34), 依次算出 , 然后把它们代入式(7.33), 算出节点xi+1上的yi+1 和zi+1值。 对于具有三个或三个以上方程的方程组的初值问题,也可用类似方法处理,只是更复杂一些而已。此外,多步方法也同样可以应用于求解方程组初值问题。 例7.6 用改进的Euler法求解初值问题 取步长h=0.1,保留六位小数。 解: 改进的Euler法公式为预报: 校正: 将 及h=0.1代入上式,得 null由初值 ,计算得 null7.7.2 高阶方程组 高阶微分方程(或方程组)的初值问题,原则上都 可以归结为一阶方程组来求解。例如,有二阶微分方 程的初值问题 (7.35) 在引入新的变量 后,即化为一阶方程组初值问题:(7.36) 式(7.36)为一个一阶方程组的初值问题,对此可用7.7.1中介绍的方法来求解。例如应用四阶龙格-库塔公式得 null(7.37) (7.38) null消去 ,上式简化为: (7.39) (7.40) 上述方法同样可以用来处理三阶或更高阶的微分方程(或方程组)的初值问题 null例7.7 求解下列二阶微分方程的初值问题 取步长h=0.1 解:先作变换:令 ,代入上式,得一阶方程组 用四阶龙格-库塔方法求解,按式(7.37)及(7.38)进行计算: 取步长 , , ,时 nullnull然后计算 时的 y2和z2;依此类推,直到i=9时的y10和z10,即可得到其数值解。 null本章小结 本章介绍了常微分方程初值问题的基本数值解法。包括单步法和多步法。单步法主要有欧拉法、改进欧拉法和龙格—库塔方法。多步法是亚当姆斯法。它们都是基于把一个连续的定解问题离散化为一个差分方程来求解,是一种步进式的方法。用多步法求常微分方程的数值解可获得较高的精度。 实际应用时,选择合适的算法有一定的难度,既要考虑算法的简易性和计算量,又要考虑截断误差和收敛性、稳定性。 null 龙格-库塔法较为常用,适用于多步方法中作初值计算和函数f(x,y)较为简单的场合。四阶标准龙格—库塔法精度高,程序简单,易于改变步长,比较稳定,也是一个常用的方法,但计算量较大。当函数f(x,y)较为复杂,可用显式亚当姆斯方法或亚当姆斯预测—校正方法,不仅计算量较小,稳定性也比较好,但不易改变步长。 一般采用龙格—库塔法提供初值y1, y2, y3,然后用亚当姆斯外推公式求得预测值 ,再由亚当姆斯内插值求得校正值yi+1,如此求得的值近似程度好且节省计算量,是一种较好的方法。null作业 习题一7.1~7.9
本文档为【数值分析-第七章常微分方程数值解】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑, 图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
该文档来自用户分享,如有侵权行为请发邮件ishare@vip.sina.com联系网站客服,我们会及时删除。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。
下载需要: 免费 已有0 人下载
最新资料
资料动态
专题动态
is_459107
暂无简介~
格式:ppt
大小:819KB
软件:PowerPoint
页数:0
分类:工学
上传时间:2012-12-09
浏览量:82