《第14章计算机数学基础下常微分方程的数值解法.docx》由会员分享，可在线阅读，更多相关《第14章计算机数学基础下常微分方程的数值解法.docx（8页珍藏版）》请在三一文库上搜索。
1、第 14 章常微分方程的数值解法常微分方程的数值解法的基本思想：求解如下形式的一阶常微分方程：dydy1f ( x, y)xdxdxy( x0 ) y0y(0)2只有少数具有较简单形式的微分方程才能求出其精确的解析解 （解析表达式表示的解）。在实际问题中往往不能求出常微分方程精确的解析解，因此只能利用数值方法 来求出微分方程的解在某些点上的 近似值，通常称为数值解法。在上述初值问题解的存在区间 a,b 内，求它在一系列节点：a=x0<x1<x2<x3 <xn=b上的近似值yk，即 yky(xk)(k=1,2,3, ,n), h= xk+1-xk(k=0,1,2, ,n-
2、1) 称为步长，一般取 h 为常数 ( 取等距节点 ) 。一、 欧拉法和改进欧拉法1欧拉公式设一阶常微分方程的初值问题：dyf ( x, y)dxy( x0 ) y0求微分方程根的几何意义：寻找一条经过点（x0,y 0）的平面曲线 y=y(x)，这条曲线在点 (x,y) 处的切线斜率是 f(x,y)。欧拉折线法的几何意义：用一系列折线来逼近所求的曲线。yx作法：过点 (x 0,y 0 ) 以 f(x 0,y 0) 为斜率作切线，其方程为：y y0 f ( x0 , y0 )( x x0 )当 x=x1 时，切线的值记为 y1, 以 y1 作为 y(x 1 ) 的近似值 , 即 y(x 1) y
3、1, 再过点(x 1,y 1 ) 以 f(x 1,y 1) 为斜率作切线，其方程为：yy1f ( x1 , y1 )( xx1 )重复以上做法，可得到切线的公式为：yynf (xn , yn )( xxn )当 x=xn+1 时，得到 y(x n+1) 的一个近似值y( xn 1 )yn 1ynf ( xn , yn )( xn 1xn )当所取的节点 xk 为等距时，有 : 步长为 h=xk+1-x ky( xk 1 )y k 1y k hf ( xk , y k )xk 1x0(k1) h欧拉公式：y( xk)ykykhf ( xk , yk )(k, , ,., n )xkxkh使用欧拉
4、公式求解一阶常微分方程初值问题近似解的方法称为欧拉法 （或欧拉折线法）。解题步骤：先明确所给的条件中的已知点(x 0,y 0 ) ，步长 h，以及所给出的f(x,y)，再套用欧拉公式中的递推公式，依次求出y1,y 2,y 3。例 1：取步长 h=0.1,用欧拉法求解初值问题yyy 的计算公式xy( )解：已知 x0 =1,y 1=1; h=0.1;x k=1+0.1k欧拉公式y(x)yyhf (x , y )xkk 1k 1kkk (k 0,1,2,.,n 1)x0kh此处 f (x, y)yy ，迭代公式为：xy k 1yky kyk )y k (1.10.1),0.1(20.1k)2(1
5、0.1k )(1(k0,1,2,., n 1), y01y( x0 ), x01例 2：用欧拉法求初值问题dy2 xdxyyy ( 0 )1在 x=0.1,0.2,0.3, ,1.0 处解的近似值。解：已知 f(x,y)= y2x , x 0=0 , y 0=1 , 取等步长 h=0.1, n=10y由欧拉公式得y(xk 1) yk 1 ykhf (xk , yk ) =yk+0.1×（ y k2 x k ）y kxk=x0+kh=0.1k（ k=0,1,2, ,9）依次求出解在各点的近似值，把计算过程列表（数值计算表） 。欧拉法的局部截断误差： 对应于节点 xk+1 处的精确解为
6、y( xk+1), yk+1 是用欧拉公式得到的近似解，则称y(xk+1) yk+1 为欧拉公式的局部截断误差。欧拉公式的22局部截断误差是O(h )，表示与 h 是同阶无穷小量，即称欧拉法具有一阶精度。例 ：用欧拉法解初值问题yyxy (x . )，取步长 h=0.2.3y()解： h=0.2,f ( xk , yk ) =yxy2 ，n=3,x k=0.2k;首先建立欧拉迭代格式：yk 1yk hf (xk , yk ) ykhykhxk yk20.2yk (4 xk yk )(k0,1,2)当 k=0，x1=0.2 时，已知 x0=0, y0=1，有y(0.2)y1=0.2 ×
7、1(4 0×1) 0.8当 k1，x2=0.4 时，已知 x1=0.2, y1 =0.8 ，有y(0.4)y2=0.2 ×0.8 ×(4 0.2 ×0.8) 0.614 4当 k=2，x3=0.6 时，已知 x2=0.4, y2=0.6144 ，有y(0.6)y3=0.2 ×0.6144 ×(4 0.4 ×0.4613) 0.82改进欧拉公式在初值问题中，对方程y'= f(x,y)在区间 x k,x k+1 上直接积分，得到 ：y( xk 1 )y( xk )xk 1f ( x, y(x) dxxk对上式可用积分梯形
8、公式，即得：y( xk 1 )y( xk )hk f ( xk , y( xk ) f ( xk 1 , y( xk 1 ) 2用 yk+1,y k 替代 y(x k+1),y(x k ) 得到计算公式：y( xk 1 ) yk 1 ykhk f (xk , yk ) f ( xk 1, yk 1)2hkxk 1 xk上式称为求 y(x k+1) 的近似值的 梯形公式。梯形公式并不能直接求得yk+1, 称为隐式形式，需要解方程才能得到。 因此一般先用欧拉公式求出初始近似值， 然后再按迭代方式求解。取等距节点，即 h 为常数，有预报值yk 1ykhf ( xk , yk )yk 1ykh也可表示
9、成下式：校正值 f ( xk , yk )f ( xk 1 , yk 1)2yk 1ykh f ( xk , yk ) f (xk 1, yk hf (xk , yk )2或表示成平均形式：yyypy khf ( x k , y k )cy khf ( x k 1 , y p )k11( y py c )2上式称为改进的欧拉公式。 用改进欧拉公式去求解一价常微分方程初值问题的近似解的方法称为改进欧拉法。改进欧拉公式的局部截断误差：y(xk 1)yk 1o(h3 )3即具有二阶精度，是h 的同阶无穷小量。例 4：用欧拉预报校正公式求解初值问题y yy sin x，取步长y( )h=0.2, 计算
10、 y(0.2), y(0.4)的近似值解：步长 h=0.2,此时 f(x,y)= y y2sin x欧拉预报校正公式为：预报值yk 1ykhf (xk , yk )校正值yk 1ykh f (xk , yk ) f (xk1 , yk1 )有迭代格式：2预报值yk 1ykh(yk yk2 sin xk )yk (0.80.2yk sin xk )校正值yk 1ykh (yk yk2 sin xk )( yk 1yk12 sin xk 1 )22yk (0.90.1yk sin xk ) 0.1( yk1yk 1sin xk 1)当 k=0，x0 =1, y0=1 时， x1=1.2 ，有yy(
11、 . ysin x)( .sin).y( .)y( .sin ). ( .sin. ).当 k，x1=1.2,y1=0.71549时， x2，有=1=1.4y 2y1 (0.8 0.2 y1 sin x1 )0.71549(0.8020.71549 sin1.2)0.47697y( . )y.( .sin . ). (.sin. )0.52608例 5：用改进的欧拉法平均公式，取步长h=0.1 ，求解初值问题y x y(0 x 0.2) y(0) 1计算过程保留 4 位小数。解 首先建立迭代格式 :x 0=0,y 0=1,x 1=0.1y pykhf ( xk , yk ) hxky k (1
12、 h)ycy khf (xk 1, y p ) yk (1 h h 2 ) hxk 1h2 xkyk1 y p yc 1h)xk hxk 1 (1h 2) yk1h(1h222当 k=0 时， x0 =0, y0=1, x1=0.1, 有y110.12 0.1 (1 0.1) 0 0.1 0.1 (1 0.1) 1 1.1122当 k=1 时， x1=0.1,y1=1.11,x2=0.2, 有y 2 1 0.1(1 0.1)0.1 0.1 0.2(1 0.10.12) 1.11 1.242 122所求 y(0.1)1.11;y(0.2)1.242 1例 6：对初值问题 yy, y( )，hn证
13、明：用梯形公式求得的近似解为ynh证明 解初值问题的梯形公式为ynyh f (xn, yn)f ( xn, yn)nf ( x, y)y整理成显式 :反复迭代，得到ynh ynhyynynynh yn yn ynhynhhhn.hynynyhhhnhh三、龙格库塔法dyf (x, y)求解此类初值问题时dx，由微分中值定理可知，存在有y(x0 )y00<<1，使得： y(xk 1)y(xk )hy ( xkh)y( xk 1 )y( xk )hf (xkh, y( xkh)f ( xkh, y ( xkh) 称为区间 xk , xk 1 上的平均斜率关 键即为寻找一种计算平均斜率的
14、计算方法，就可求得近似值yk+1。若仅取 xk 处的斜率 f(x k ,y k) 作为平均斜率的近似值，可得：y(xk 1)y(xk )hf (xkh, y(xkh)yk 1ykhf ( xk , yk ) （欧拉公式）若取 xk 处的斜率 f(x k ,y k) 和 xk+1 处的斜率 f(x k+1 ,y k+1) 的平均值，作为平均斜率的近似值，可得：yk 1ykh f (xk, y )f ( xk 1, yk 1)2k（梯形公式，进一步可得改进欧拉公式）一般取 N 个点的斜率的某种线性组合作为平均斜率：Nf ( x kh , y ( x kh )m mm 1其中 m表示第 m个点的斜率
15、值的大小，m表示该点斜率所对应的 系数 。则计算公式变为：Nyk 1yk hm mm 1当 N 取为 2 时，所得公式称为两点的二阶龙格库塔公式；当 N 取为 3 时，所得公式称为三点的三阶龙格库塔公式；当 N 取为 4 时，所得公式称为四点的四阶龙格库塔公式；显然当取的点越多时所得的近似值 yk+1 越接近精确值，即其局部截断误差变小、精度越来越高，但相应的计算量也越来越大。龙格库塔法的局部截断误差：3二阶龙格库塔法为O(h )，具有二阶精度；4三阶龙格库塔法为O(h )，具有三阶精度；5四阶龙格库塔法为O(h )，具有四阶精度；四阶龙格库塔法公式：hyk 1y k(12 22 34 )6其
16、中 h 为步长，且有：1=f(xk,yk)；2=f(xk+ 1 h，yk+h1)；23=f(xk+ 1 h，yk+h2)； 4=f(xk+h，yk+h 3)2四阶龙格库塔法的优点：1>它是一步法； 2.>精确度高。缺点：计算量大。yy例 7：写出用四阶龙格库塔法求解初值问题y ( )的计算公式，取步长 h=0.2，计算 y(0.4)的近似值，计算过程保留4 位小数。解：此处 f(x,y)=83y, 四阶龙格库塔法公式为hyk 1yk6( 122 234 )其中1=f(xk,yk)；1h，yk+h 1)；2=f(xk+23=f(xk+1h， yk+h2)； 4=f(xk+h， yk+h 3)2本例计算公式为：ykyk. ()其中1=83 yk；2=5.62.1 yk；3=6.322.37yk;4=4.2081.578ykyk 10.2(83yk2(5.62.1yk )2(6.32 2.37yk )(4.208 1.578yk )yk61.20160.5494yk (k0,1,2,.,n1)当 x0=0,y0=2,y( . )y.y.y( . )y.y.