0101 第五章 常微分方程的数值解法

第五章常微分方程的数值解法许多工程实际问题的数学模型可以用常微分方程来描述，但是，除了常系数线性微分方程和一些特殊的微分方程可以用解析方法求解以外，在生产实践和科学研究中遇到的绝大多数微分方程往往比较复杂，不能给出解的解析表达式；有时候即使能用解析表达式来表示，又因为计算量太大而不实用，有些是已经有了求解的基本方法的典型方程，但实际使用时也是有困难的。很多实际问题只需要求某些点上的函数值，不必求出函数表达式。因此，我们有必要研究常微分方程的数值近似解法。本章主要研究一阶常微分方程初值问题的几个常用解法。例：某伞兵与降落伞质量为90kg，当伞张开时以12m/s速度垂直下降。假设空气阻力与伞的下降速度成正比，在速度为6m/s时，测得的空气阻力为353N。试求伞兵开伞后，第1秒末、第2秒末、直到第5秒末时各时刻的速度。解：设阻力为，则fvf8586353.由力学知识可得，vgdtdv8589090.120)(v50t289smg/.vv158989..120)(v50t本章主要内容一阶常微分方程初值问题欧拉方法、梯形格式、改进的欧拉方法龙格－库塔（Runge-Kutta）方法第一节引言一阶方程的初值问题：0yaybaxyxfdxdyy)(],[),(由常微分方程理论，在满足一定条件时存在唯一解函数。),(yxf)(xy要计算出解函数y(x)在一系列节点a=x0x1…xn=b处的近似值。),...,()(nixyyii0常微分方程的数值解节点间距为步长，通常采用等距节点，),...,(101nixxhiii即取hi=h(常数)。在这些节点上采用离散化方法，（通常用数值积分、微分、泰勒展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解作为的近似值。这样求得的就是上述初值问题在节点上的数值解。一般说来，不同的离散化导致不同的方法。)(ixyiyixiy第二节欧拉方法1、欧拉方法的导出0yaybaxyxfdxdyy)(],[),(对于上面的初值问题，若将函数在处的导数用差商来表示，则有，)(xyix)1,...,0(),(1niyxfhyyiiii),(0001yxfhyy),(1112yxfhyy),(111nnnnyxfhyy…xxy21)(解：根据欧拉公式可以得到：此外，方程的真解为：2、举例例1用欧拉法求初值问题取步长h=0.2。1)0(2yyxydxdyy10201)()(*),(yyyxyhyyxhfyyiiiiiiiii012345xiyi01.00000.20001.20000.40001.37330.60001.53150.80001.68111.00001.8269求解结果如下：y(xi)1.00001.18321.34161.48321.61251.7321i=y(xi)-yi00.01680.03170.04830.06860.0949()yyxxy0x1x2x0P1P2P根据已知条件：曲线y(x)上的点（x0,y0）及该点处曲线的导数f(x0,y0)，则可以得到过该点的直线：))(,(0000xxyxfyy该直线与x＝x1的交点P1，则P1的纵坐标y1为：))(,(010001xxyxfyy),(000yxhfy就用y1作为y(x1)的近似值…逐次进行后可以得到一条折线P0P1…Pn，该折线看作是初值问题的积分曲线的近似，因此欧拉方法也称为欧拉折线法。3、欧拉方法的几何意义欧拉法虽然形式简单，计算方便。但从上述几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙,精度也低，即误差太大。特别当的曲线曲率较大时，欧拉法的效果更差。二、梯形格式11iixxiidxxyxfxyxy))(,()()(对微分方程两边积分得，))(,())(,(iixxxyxhfdxxyxfii1)1,...,0(),(1niyxfhyyiiii欧拉方法左矩形公式))(,())(,(111iixxxyxhfdxxyxfii11iixxiidxxyxfxyxy))(,()()(),...,(),(10111niyxfhyyiiii后退的欧拉方法右矩形公式)))(,())(,(())(,(1121iiiixxxyxfxyxfhdxxyxfii11iixxiidxxyxfxyxy))(,()()(),(),(1112iiiiiiyxfyxfhyy梯形方法梯形公式和欧拉公式相比较，梯形公式精度有所提高；但梯形公式是隐式的，在计算yi+1时候也只用到前一步的值yi，若yi已知，将yi代入公式求解时，一般不能直接得到yi+1，而需要通过其他方法（比如迭代法）求解，计算量较大。实际中，将欧拉公式和梯形公式联合使用，从而得到改进的欧拉方法。三、改进的欧拉方法将欧拉公式和梯形公式联合使用，先用欧拉公式得出一个y(xi+1)的粗糙近似值，称为预估值，然后对预估值使用梯形公式对它进行精确化，得到较为精确的近似值yi+1，称之为校正值，计算公式为：),(iiiiyxhfyy1),(),(1112iiiiiiyxfyxfhyy这样的预估校正系统称为改进的欧拉方法112()2iihyykk),(iiyxfk1),(iiiiyxhfyy1),(112iiyxfk为了便于编写程序，常将上面的公式改写为如下式：改进的欧拉方法与梯形方法具有同样的精度级别，但改进的欧拉方法为显式格式，使用更方便。xxy21)(解：根据改进的欧拉公式得：此外，方程的真解为：举例例1用改进的欧拉法求初值问题取步长h=0.2。1)0(2yyxydxdyyiiiiiyxyhyy21),(),(1112iiiiiiyxfyxfhyy11222iiiiiiiyxyyxyhyi01234501.00000.20001.18670.40001.34830.60001.49370.80001.62791.00001.7542xiyi求解结果如下：y(xi)1.00001.18321.34161.48321.61251.7321i=y(xi)-yi00.00350.00670.01050.01540.0222例：某伞兵与降落伞质量为90kg，当伞张开时以12m/s速度垂直下降。假设空气阻力与伞的下降速度成正比，在速度为6m/s时，测得的空气阻力为353N。试求伞兵开伞后，第1秒末、第2秒末、直到第5秒末时各时刻的速度。vv158989..120)(v50t012.00001.000013.96002.000014.63953.000014.87504.000014.95675.000014.9850012.00001.000013.31972.000014.05893.000014.47294.000014.70485.000014.8346欧拉法改进的欧拉法第三节龙格－库塔方法一、龙格－库塔法的基本思想拉格朗日中值定理：在闭区间[a,b]上连续，开区间(a,b)上可导，则至少存在(a,b)上的一点，使得下式成立：)(xfabafbff)()()(1、预备知识2、推导)()()(yxxxyxyiiii11)()()(yhxyxyii1),(1iixx)()()(yhxyxyii1Khxyxyii)()(1))(,()(yfyK称为在上的平均斜率。)(xf],[1iixx问题转化为如何对进行数值计算。K特例：1）取一个点上的斜率值作为平均斜率的近似值。ix),(iiyxf*K),(iiiiyxfhyy1欧拉公式精度较低2）用两个点和上的斜率和的算术平均值作为平均斜率的近似值。ix1ix1K2K*K),(),()(11212112hKyxfKyxfKKKhyyiiiiii改进的欧拉公式精度稍高根据前面的分析得到：如果在区间上多预测几个点的斜率值，然后取其加权平均作为平均斜率的近似值，有可能构造出具有更高精度的计算公式，此即为龙格库塔算法的基本思想。],[1iixx二、龙格－库塔公式),(),()(12122111phKyxfKyxfKKKhyyipiiiii1、二阶龙格－库塔公式选择参数使得算法具有2阶精度p,,21211221p满足该式的参数不止一组，而是一簇，所有满足条件的公式通称为二阶龙格库塔公式。改进的欧拉公式是二阶龙格库塔公式的其中一种。21121＝＝,p2、四阶龙格－库塔公式二阶的龙格－库塔公式是由两个点的斜率加权平均（或线性组合）得到的，同样，用四个点上的斜率线性组合就可得到四阶的龙格－库塔公式。经典的四阶的龙格－库塔公式：),(),(),(),()(3422231222143216122hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii举例例1用经典龙格—库塔公式求初值问题取步长h=0.2。1)0(2yyxydxdyy),(),(),(),()(3422231222143216122hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii求解结果如下：01.00000.20001.18320.40001.34170.60001.48330.80001.61251.00001.7321i012345xiyiy(xi)1.00001.18321.34161.48321.61251.7321i=y(xi)-yi00.00000.00010.00010.00000.0000四阶龙格—库塔法的精确度高于欧拉法和改进的欧拉法。必须说明，由于四阶龙格—库塔法有较高的精确度，大家自然会想到推导更高阶的龙格—库塔公式。但实践证明，高于四阶的龙格—库塔公式，不但计算量大，而且精度并不一定会提高，有时甚至还会降低。本章内容小结1、了解数值求解微分方程的基本概念和基本思想；2、熟练掌握一阶常微分方程初值问题数值求解的基本方法欧拉方法，改进的欧拉方法3、了解龙格－库塔方法的基本思想，并熟练掌握经典4阶龙格－库塔方法