推荐精品级实习指导书.doc

推荐精品级实习指导书 就是建模和解模两件事！建模的基础是微积分知识。 而由于微分方程求解析解的困难，基于计算机的数值计算方法就成为解模的主要手段。 而微积分知识又归结为求导公式和求积公式这两个公式。 这是此次大会于我的最大收获！正切合了此次常微分方程课程实习的指导思想！在常微分方程课程中，我们学习了几类可求解的常微分方程（组）的解法及其相关理论基础。 但是我们更知道，更多的常微分方程是不能用初等方法求解的。 然而在科学技术、工程应用等领域，常微分方程是经常遇到的主要数学模型。 一般来说，找出解的解析表达式极其困难，对大部分问题甚至是不可能的。 因此，寻求近似解法就非常必要，而微分方程数值解法是目前普遍采用的微分方程近似解法。 本实习指导书即给出常微分方程数值解法利用计算机依托数学软件（如MATLAB）求常微分方程的近似解，即给出理论解在一些离散点上的近似值。 本实习指导书中讨论的是常微分方程的初值问题00(,)f x y()dydxy xy?)2()1(如无特别说明，总认为这个初值问题的解是适定的，即解存在、唯一且连续依赖于初始条件。 要求理解求微分方程初值问题的数值解方法，包括欧拉方法和改进的欧拉方法及龙格-库塔方法。 能应用欧拉方法和龙格-库塔方法求微分方程的初值问题的数值解并画出图形。 常微分方程课程实习指导书22欧拉（Euler方法）2.1欧拉格式数值求解常微分方程的初值问题的最简单的方法是（向前）欧拉法（Euler）。 该方法的理论推导如下由),(yxfdxdy?，两边在,00hxx?上同时积分得0000()()y x(,()f xy xdxxhxy xh?而求?hxxdxyxf00),(的近似值的最直接的方法是用),(00yxf代替),(yxf（即定积分近似计算中的用矩形近似曲边梯形），于是有?),()()(0000yxhfxyhxy?用0y表示)(0xy，1y近似代替hxxhxy?010),(，则),(0001yxhfyy?依此令101,?nnynhxx近似代替)(1?n xy，则?)(2,1,0),(001xyynyxhfyynnnn? （3）?1称为步长。 这就是欧拉法的计算公式，nnxxh?一般而言，并不要求步长相等，则有?),()(11nnnnnnyxfxxyy （4）例1以1.0?h为步长，用欧拉法求初值问题?1)0(yyxedxdyx的数值解，并与精确解xexxy)2 (21)(2?比较。 解由欧拉法有?y?12,1,0)(01nyexhyynxnnnn?在使用欧拉法数值的求解过程中，我们发现计算过程非常简单，即由0y可直接计算1y，由1y可直接计算?,2y，无需用迭代方法求解任何方程，因此也称为显式格式。 一般而言，我们难以得到精确解)(n xy，用欧拉法数值求得的近似解n y，误差为|)n(|nxyy?；另常微分方程课程实习指导书3一方面，不论何种格式都不能利用计算机求得它们的精确解n y，这是因为计算机运算，总是有限位二进制运算，因此对十进制的数学运算总会出现舍入误差及其在计算过程的传递。 因此计算机输出的是欧拉法的近似解()(nnnnyyxyy?n y，而不是精确解n y。 由于)()nnxyy? （5）可见，为了使计算得出的解n y是)(n xy好的近似。 我们要求欧拉法的精确解n y是微分方程精确解)(n xy的很好近似，特别要求当步长h充分小时，所得的近似解y?n y是n y能足够精确地逼近精确解)(n xy。 换言之，要求0?h时，).(nnxyn y好的近似。 由于计算机在计算过程会不断地产生舍入误差，本问题的讨论相当复杂。 为了简化讨论，我们设想计算机对欧拉格式计算过程完全精确，每步都没有误差，|nnyy?的值完全由决定。 要求因此|y?|00yn y是n y的好的近似则相当于要求欧拉格式解对初始值具有连续依赖性，这种解对初值的连续依赖性就称为稳定性。 问题 （1）称为格式的收敛性问题。 问题 （2）称为格式的稳定性问题。 格式的收敛性、稳定性研究微分方程数值解法最基本的理论研究工作，具有重要的实用意义。 一个算法格式只有当它是稳定的、收敛的，才是能用的、有意义的。 2.2收敛性研究简介前面已指出，收敛性问题，即研究?nnxy?)(,00xy，yxnhxhn?时的问题。 记111)(?ny为欧拉方法的整体截断误差，而由)(n xy利用欧拉公式得到的1?ny与)(1?n xy之差?1ne|)1?(|1?nnxyy称为局部截断误差，也可记为nR。 定理1假定,)(02XxCxyy?，则欧拉方法的局部截断误差nR满足221|MhRn?其中h为步长，),(yx关于y满足Lipschitz条件，L为相应的Lipschitz常数，则欧拉|)x(|max0xx?yMX?。 定理2设f方法的整体截断误差n?满足)1?(|) (0)(00?LxXLxXnehLRe?其中R为局部截断误差的上界。 ,(xf定理3设)y关于y满足Lipschitz条件，L为相应的Lipschitz常数，,)(02XxCxyy?且当0?h时，)(00xyy?，则欧拉方法的解n y一致收敛到初值问题常微分方程课程实习指导书4的解)(n xy，并有估计式整体截断误差n?满足)1?(|) (0)(00?LxXLxXnehLRe?其中R为局部截断误差的上界。 )1?(2|) (0)(00?LxXLxXneLMhe?如果)(00xyy?，即欧拉方法的整体截断误差n?与h同阶，由nR的表达式可知，)(2hoRn?，这说明局部截断误差比整体截断误差高一阶。 我们称欧拉方法一阶格式。 2.3稳定性的研究简介稳定性问题是决定欧拉法在利用计算机计算中能否得到精确解的关键性问题，只有稳定的算法才可能是有用的算法。 定义1如果存在正整数c及0h，使对任意初始值,00zy用?y?01),(yxhfyynnnn)n与?z?01,(zxhfzznnn计算所得之解?nnzy,满足估计式?0000;0|,|xXnhhhzyczynn?则称欧拉法是稳定的。 注意，这里分别是以00,zy为初值得到的精确值，毫无舍入误差，因此这里的稳定性定义是对初值的稳定性，即研究初值误差在计算过程的传递问题。 定理4在定理2的条件下，欧拉方法是稳定的。 由定理2，我们看到初始误差00?，则整体截断误差的阶完全由局部截断误差的阶所决定。 因此要提高数值算法的精度，往往从提高局部截断误差的阶入手，这也是构造高精度数值方法的主要依据。 常微分方程课程实习指导书53梯形法、隐式格式的迭代运算用梯形公式近似计算积分，则)(,()(,(2)(),()()(11111?nnnnxxnnnnxyxfxyxfxxdxyxfxyxynn因此有)(,()(,() (21)()(1111?nnnnnnnnxyxfxyxfxxxyxy可得梯形公式为),(),()(211111?nnnnnnnnyxfyxfxxyy显然这是一个隐式格式。 它的局部截断误差阶为)(3ho，较之欧拉法高一阶。 类似于欧拉法，对梯形法可以平行地建立它的整体截断误差阶为)(2ho，以及格式的收敛性稳定性等定理。 梯形法是一个隐式格式),(),(21111?nnnnnnyxfyxfhyy如何求解1?n y呢？我们采用迭代法，其格式如下 (1)()1?11 (0)1?12(h f xy,)(,)ppnnnnnnnyyf xyy?初始猜测易证梯形法迭代格式收敛的充分条件是12?hL。 在实际计算中，令0?p有下面的预报-校正格式?1校正格式预报格式),(),(2),()0 (111)0(1nnnnnnnnnnyxfyxfhyyyxhfyy当然也可以迭代多次?校正格式预报格式),(),(21),() (11)1 (1)0(1pnnnnnpnnnnnyxfyxfhyyyxhfyy当步长取得适当小，用预报格式已能算出比较好的近似值，故迭代收敛很快，通常只需迭代两三次就可以满足精度要求，如果迭代多次仍不收敛，说明步长过大，必须减小步长，再进行计算。 梯形法较之欧拉法提高了精度，但增加了迭代次数，因此增加了计算工作量。 例2用预报-校正格式解初值问题常微分方程课程实习指导书6?|?11,010xyxxyy取1.0?h。 解预报-校正格式为?)1()1 (21)1 (1)0 (11)0(1nnnnnnnnnnxyxyhyyxyhyy也可以写成?与欧拉格式的结果比较，可以看出，欧拉法精确度较低，预报-校正格式精确度有所?9,3,2,1,0)1(1.0)1(1.0)(21112121x1?nxkykykkkyynnnnnn改善。 可选择的问题一显式Euler方法是解常微分方程初值问题最直观、最简单的方法。 考虑一个简单的问题()()sin (0)?0y tytty?由于该方程是线性的，故它的解析解易于得到。 1.Matlab实现Euler方法。 在0,10t?的范围内用Euler方法求解上述问题，并将求出的解用matlab的图形工具画出来，同时将方程的解析解也画在同一张图上。 2.选择不同的步长h（至少取3个有区分度的值），比较不同的步长所得数值解逼近原微分方程解的精确程度。 请各自选择其他的方程尝试上述算法，要求不同同学的方程不能一样。 自变量的范围也可调整，要求较多的节点（才能看出区别，才能说明问题）。 问题二考虑VanderPol方程23213()0d udtdudtdudtu?。 该方程可以写成如下的一阶方程组的形式313()dxdtuxxdudtx?经过变量替换/u?得到313()dxdtxpxddtx?常微分方程课程实习指导书7其中p为引入的模型参数，当p-1时即前方程。 这是一个最简单的心脏模型，当参数p取不同值时，解的波形也在变化。 写一个简单的matlab程序求解它。 4一般单步法、Runge-Kutta格式前面，我们研究了欧拉法和梯形法，它们有一个共同的特点，即在格式中只包括,?nnnyxy的值，或者说由，仅使用11,nx1?nnxxn y的值计算出1?n y的值，这种格式称为单步格式，下面研究一般单步法。 4.1一种构造单步法的方法泰勒级数法初值问题00(,)fxy()dydxy xy?的解)1?)(qxy阶可微，将)(xy在0x处展开成泰勒级数，有)()(!)(!2)()() (10)(02000?qqqhoxyqhxyhxyhxyhxy?由方程可得xy?)(yyxfyxfyxfxyy)xfy?x?),(),(,(),()(000因此(y?)() (2)()(),(),(),()20000000yxy?yy?xy?xx?y?xfffffffffxydxdxyyxfyxfyxfx?式中，?,xx?y?xfff?都是),(yxf相对于变量的偏导数。 于是有)(),(,()()(10000?qhohxyxhxyhxy?其中)(,(11110000|)(,(!1j(),(,(xyxjjjqjxyxfdxdhhxyx?舍去)(1?qho，可得),(0001hyxyy?),(1112hyxyy?yn?),(1hyxynnn?h称),(1yxyynnnn?为一般单步法，显然局部截断误差)()(1?qnnnhoyxye常微分方程课程实习指导书84.2一般单步法基本理论推广由泰勒级数法得到的q阶单步法，下面我们给出一般单步法基本理论。 定义2给出单步法?),(1hyxyynnnn?),(,(hxyx为任意关于),(,(hxyx的函数，其对于微分方程),(yxfdxdy?的解)(xy满足)(),(,()()(1?qhohxyxhxyhxy?且q为使上式成立的最大整数，则称?yxyynnnn?,2,1,0),(1?nh为q阶单步法，由)(),(,()()(1?qhohxyxhxyhxy?得)0,y,()()(limh?0xhxyhxy?因此有如下定义。 定义3如果),()0,y,(yxfx?，则称单步法),(1hyxyynnnn?为与初值问题相容的。 定理5如果),(hyx?对于000,hhXxx?以及所有实数y满足Lipschitz条件，则单步法稳定。 证明略下面建立稳定性与收敛性的关系。 前面已建立了单步法稳定性定理，现在的问题是其精确解当0?h时，n y是否收敛到)(n xy，我们利用格式的相容性把两者联系起来。 ?以及所有实数关于定理6如果),(hyx?对于?000,hhXxx?hyx,满足Lipschitz条件，则)0,(fyx?),(1hyxyynnnn?收敛的充要条件是格式相容，即满足),(yx。 证明略关于单步法的整体截断误差n?，我们有定理7在定理5的条件下，如果局部截断误差?的整体截断误差nR为)(1?qho，则单步法),(1hyxyynnnnn?满足)1?(|) (0)(00?LxXqLxXnehL?ce?特别地，若00?有)(qnho?，整体截断误差比局部截断误差低一阶。 证明略。 常微分方程课程实习指导书94.3龙格库塔格式从前面讨论可见，构造单步法的关键在于构造),(hyx?，使)(),(,()()(11?qnnnnhohxyxhxyxy?中的局部截断误差阶尽可能提高。 前面我们利用泰勒级数法构造了一个欧拉方法，这时,1),()(,()()(21?qhoxyxhfxyxynnnn局部截断误差与)(2ho同阶，这是一个一阶格式。 为了要求2?q，利用泰勒级数法得到一个二阶格式),(),(),(2),(21nnnnynnxnnnnyxfyxfyxfhyxhfyy?这时我们有)()(,()(,()(,(2)(,()()(321hoxyxfxyxfxyxfhxyxhfxyxynnnnynnxnnnn?格式计算过程中要求函数的两个偏导数在)(,(nnxyx处的值，比较麻烦，可以预计，利用泰勒级数法得到的高阶格式需要求更多的偏导数值，计算繁复。 那么是否可以避免计算偏导数，),(hyx?呢？分析梯形法的预报校正格式而得到高阶单步格式的?1校正格式预报格式),(),(2),()0 (111)0(1nnnnnnnnnnyxfyxfhyyyxhfyy这时可以写成?),(1hyxhyynnnn?其中),(, (21), (21),(11nnnnnnnyxhfxfyxfhyx?记),(,2111nnyxfKc?),(,21122hKyhxfKn?则21212x1),(KKhyxnn?这时从1?nnx，单步法可以分两级进行；K?第一级计算),(1nnyxf第二级计算),(12hKyhxfKnn?最后计算)(22111KcKchyynn?。 常微分方程课程实习指导书10可以证明这个格式的局部截断误差阶为)(3ho，我们称之为二级二阶的龙格库塔方法。 一般的二级二阶的龙格库塔格式可以写成?),(),()(121221221y11KhbyhaxfKxfKKcKchyynnnnnn适当地选择参数21221,ba，使局部截断误差)(31hoRn?因此要求满足?这是一个含四个参数、三个方程的方程组，因此有一个自由参数，解答不唯一。 ?212112212221cbca （1）1,21,2121221?ba，即得二级二阶的龙格库塔方法。 ?x?),(),y()(21121211hKhfKyxfKKKhyynnnnnn1 （2）2,1,021221?ba，即得?)2,23(),(12121hKyhxfKyxfKhKyynnnnnn2 （3）3,4,4121221?ba，即得?类似可得三级三阶的龙格库塔的一般算式?)32,32(),()3(4f1121211hKyhxfKyxKKKhyynnnnnn?h?K?),(),()K,hb()(232131331212213322111KhbKhbyhaxfKyaxfKyxfKcKcKchyynnnnnnnn常微分方程课程实习指导书11适当地选择参数323212321,31bbabac，使局部截断误差)(41hoRn?因此要求满足?616131) (31) (3121)(21121323232332323122213231332122322222323132213222321bbcabbbcbbbcabaacabbbcacac由上方程组中最后二个方程的?，于是上方程组可化为212ba?，代入第二个方程与第三个方程联立，得32313bba?这是八个参数、六个方程的方程组，有两个自由参数，故有无穷多个解，特例如下1,66?613121121323322222322232313212321bbacacacabbabac （1）,1?,2,21,1,24,131322132221?bbbaac得三级三阶的龙格库塔算法?)2,()2,2h(),()4(62131213211hKhKyxfKKhyhxfKyxfKKKKhyynnnnnnnn （1）,0,32,31,32,31,32,0,4131322132321?bbbaac得Heum三级三阶的龙格库塔算法常微分方程课程实习指导书12?类似可得四级四阶的龙格库塔的一般算式?)32,32()3,3()h,()3(42312131y1KhyhxfKKyhxfKxfKKKhyynnnnnnnn?如上推导，为了达到四级四阶格式，可得13个参数满足11个方程?212bba?f?),(),(),(),()(343242141432321313312122144332y2x111hKbKhbKhbyhaxfKKhbKhbyhaxfKKhbyhaxfKKKcKcKcKchyynnnnnnnnnn?c?c?24112181)(61413121143322443234422243222344334224323234334422432233443333222442332223333223321434241432313bbacbacbacbacababacbaacbacbacbacacacacacacacacacacbbbaba该方程中有两个自由参数，下面给出二组解经典四级四阶龙格库塔格式取定2132?aa，则得常微分方程课程实习指导书13?),()2,2h()2hK,2h(),(h)22(63423121432y11hKyxfKyxfKKhyxfKxfKKKKKhyynnnnnnnnnn这是最为著名的经典四级四阶龙格库塔格式。 （2）取定516387,4.032?aa，则得另一四级四阶龙格库塔格式?)?)83286432.3050964.321810038.0,()hK15875966.0hK29697760.0hK,455737254.0yhn()4.0,4.0(,()17118478.02055354.155148053.0K17476028.o(321421312143211