第七章 常微分方程数值解法
常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。
第一节 欧拉法
求解常微分方程初值问题
?????==0
0)()
,(y x y y x f dx
dy
（1）
的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<< {}n y 称为问题的数值解，数值解所满足的离散方程统称为差分格式，1 --=i i i x x h 称为步长，实用中常取定步长。 显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。 定理 设函数()y x f ,在区域+∞≤≤-∞≤≤y b x a D ,: 上连续，且在区域D 内满足李普希兹(Lipschitz)条件，即存在正数L ，使得对于R 内任意两点()1,y x 与()2,y x ,恒有()()2 121,,y y L y x f y x f -≤-则初值问题(1) 的解()x y 存在并且唯一。 一、欧拉法(欧拉折线法) 若将函数 )x y 在点n x 处的导数()n x y ' 用两点式代替， 即 ()h x y x y x y n n n ) ()(1-≈ '+，再用n y 近似地代替()n x y ,则初值问题(1)变为 ?? ?==++=+ ,2,1,0),() (00 1n x y y y x hf y y n n n n (2) (2)式就是著名的欧拉(Euler)公式。以上方法称 为欧 拉法或欧拉折线法。 欧拉公式有明显的几何意义。从几何上看，求解初值问题(1)就是xy 平面上求一条通过点()00,y x 的曲线()x y y =，并使曲线上任意一点()y x ,处的切线斜率为 ()y x f ,。 欧拉公式的几何意义就是从点 ) 000,y x P 出发作一斜率为()00,y x f 的 直线交直线1x x =于点()111,y x P ，1P 点的纵坐标1y 就是()1x y 的近似值；再从点1P 作一斜率为()11,y x f 的直线2x x =交直线于点()222,y x P , 2P 点的纵坐标2y 就是的近似值()2x y ；如此继续进行，得一条折线 210P P P 。该折线就是解()x y y =的近似图形，如图7-1。 图7-1 欧拉法的其它几种解释： 1． 假设()x y 在n x 附近展开成泰勒级数 ()()()()()()()() +''+ +=+''+'+=+n n n n n n n n x y h x y x hf x y x y h x y h x y x y 2 ,2 22 1 取h 的线性部分，并用n y 作为()n x y 的近似值，得 ()n n n n y x hf y y ,1+=+ 2． 对方程 () y x f dx dy ,=两边从n x 到1+n x 积分，得 ()dx x y x f x y x y n n x x n n ? += -+1 ))(,()(1 (3) 用 矩形公式计算上式右侧积分，即 ()()x x x x x d x y x f dx x y x f n n n n ? ? ++≈ 1 1 ,))(,( 并用n y 作为的近似值()n x y ，得()n n n n y x hf y y ,1+=+ 故欧拉法也称为矩形法。 二、改进的欧拉法(梯形法) 欧拉法形式简单，低，特别当曲线y=y(x)计算方便，但精度比较的曲率较大时，欧拉法的效果更差。为了达到较高精度的计算公式，对欧拉法进行改进，用梯形公式计算(3)式 右侧积分，即 ()()()()[] 11,,2 ))(,(1 +++≈ ? +n n n n x x x y x f x y x f h dx x y x f n n 并用y n 作为y(x n )的近似值，得到改进 的欧拉公式 ()()[]111,,2 +++++ =n n n n n n y x f y x f h y y (4) 上述方法称为改进的欧拉法，也称为梯形法。 不难发现，欧拉公式 是 关于y 1+n 的显式，即只要已知y n ,经过一次计算便可得y 1+n 的值，而改进的 欧拉公式是以y 1+n 的隐式方程给出，不能直接得到y 1+n 。隐式方程(4)通常用 迭代法求解，而迭代过程的实质是逐步显式化。 先用欧拉 公式 () () n n n n y x hf y y ,01+=+ 给出y 1+n 的迭代初值，然后 再用改进的欧拉公式(4)进行迭代，即有 [] ?? ? ????=++=+=+++++ ,2,1,0),(),(2),()(11)1(1) 0(1k y x f y x f h y y y x hf y y k n n n n n k n n n n n (5) 迭代过程进行到连续两次迭代结果之差的绝对值小于给定的精度ε即 () ε k n k n y y 111+++- 为止，这时取 () 11 1+++=k n n y y 然后再转入下一步计算。 下面讨论(){}k n y 1 +是否收敛；若收敛，它的极限是否满足(4)式。 假设f(x,y)满足李普希兹条件 ()()()2121,,y y L y x f y x f -≤- 则 () () () ()() () ()() () ()() 01 112 2 1 112 11 111 111111222 ,,2 ++-+-+-++-+++++++-? ? ? ??≤≤-?? ? ??≤-≤ -= -n n k n k n k n k n k n n k n n k n k n y y hL y y hL y y hL y x f y x f h y y 由此可见，只要 12 k →∞时，有0 2→?? ? ??k hL ， 所以() {}k n y 1 +收敛。又因为f(x ,y)对y 连续，当k →∞时，对等式 () ()() ( )[] k n n n n n k n y x f y x f h y y 1 111 ,,2 ++++++ = 两端取极限， 得 ()()[] 111,,2 +++++ =n n n n n n y x f y x f h y y 因 此，只要步长h 足够小，就可保证(){}k n y 1 +收敛到满足(4)式的1 +n y 。 三、预估一校正法 改进的欧拉公式在实际计算时要进行多次 迭代，因而计算量较大。在实用上，对于改进的欧拉公式(5)只迭代一次，即先用欧拉公式算出y 1+n 的预估值y ) 0(1 +n ，再用改进的欧拉公式(4)进行一次迭代得到 校正值y 1+n ，即 ()()[]?? ? ? ?=++=+=++++ ,2,1,0,,,2),(111) 0(1n y x f y x f h y y y x hf y y n n n n n n n n n n (6) 预估—校正公式也常写成下列形式： () ,2,1,0,,),(21211 212 11 =??? ? ?? ? ++==++=+n k y h x hf k y x hf k k k y y n n n n n n (7) 四、公式的截断误差 定义 若某种微分方程数值解 公式的截断误差是O(h 1 +k )，则称这种方法是 k 阶方法。为了简化分析，在进行误差分析时，我们假设前一步的结果是准确的，即在 y n =y(x n )的前提下，考虑用y 1+n 作为y(x 1+n )的近似值而产生的截断误差，这种误差称为局部截断误差。由泰勒公式 ()()()()() +''+ '+=+=+n n n n n x y h x y b x y h x y x y ! 22 1 对于欧拉公式，有 ()()()n n n n n n x y h x y y x hf y y '+=+=+,1 于是 ()()2 11h O y x y n n =-++ 则欧拉公式的截断误差为O(2 h ),所以 欧拉法是一阶方法。 对于预估—校正公式，有 ()() ()()()()()()()()()[ ] ()()()()()()()[] +'++=+++=++=++='==n n y n n n x n n n n y n n x n n n n n n n n n x y x f x y x y x f h x y x hf x y x f k x y x hf x y x f h k x y h x hf k y h x hf k x y h y x hf k ,,,,,,,,,2 11121 而 ()()() ()()()()()() x y x f x y x y x f x y x y x f x y y x ,,,?'+=''=' 于是 ()() +''+'=n n x y h x y h k 2 2 因此 ()()() +''+ '+=+ + =+n n n n n x y h x y h x y k k y y 2 2 1212 211 所以y(x 1+n )- y 1+n = O(3 h )则预估—校正公式的截断误差为O(3 h ),也即预估— 校正法是二阶方法。 可以证明，改进的欧拉公式与预估—校正公式的截断误差相同，均为O(3 h )。这里略去证明。 例 1：在区间［0,1］上以h=0.1为步长，分别用欧拉法与预估— 校正法求初值 问题??? ??=-=1 )0(2y y x y dx dy 的数值解 解：该方程为贝努利方程，其精确解为x y 21+=欧拉公式的具体形式为 ???? ? ?-+=+n n n n n y x y h y y 21 预估—校正公式的具体形式为 ()??? ? ???? ? ???? ??++-+=-=++=+11212 112)2(2121k y h x k y hf k y x y hf k k k y y n n n n n n n n 取步长1,0,1.000===y x h ,计算结果见下表：表7-1 表7-1 近似解与准确解比较，欧拉法的结果大致只有两位有效数字，而预估—校正法的结果则 有3位有效数字。 第二节 龙格—库塔法 前面讨论的欧拉法与改进的欧拉法都是一步 法，即计算y 1+n 时，只用到前一步值。龙格—库塔(Runge-Kut ta)法(简称为 R-K 方法)是一类高精度的一步法，这类方法与泰勒级数法有着密 切的关系。 一、泰勒级数法 设有初值问题 ?????==0 0)(),(y x y y x f dx dy 由 泰勒展开式 ()()()()()() ()() 1 2 1! ! 2++++ +''+ '+=+=k n k k n n n n n h O x y k h x y h x y h x y h x y x y （1） 若令 ()()() ()n k k n n n x y k h h x y h x y y ! ! 22 1+ ++ '+=+ (2) 则 ()( )1 11+++=-k n n h O y x y 即公式(2)为k 阶方法。 从理论上讲，只要解y(x)有任意阶导数，泰勒展开方法就可以构造任意阶求 y 1+n 的公 式，但由于计算这些导数是非常复杂的。如 ()()()() yy y y x xy xx y x y x f f ff f f ff f y ff f f y f x y f x y x f x y 2 2 2,++++=''+='+=''==' 所以这种方法实际上不能用来解初值问题。 二、龙格 —库塔方法(R-K 方法) R-K 方法不是通过求导数的方法构造近似公式，而是通过计算 不同点上的函数值，并对这些函数值作线性组合，构造近似公式，再把近似公式与解的泰勒展开 式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。 我们先分析欧拉法与预估—校正法。对于欧拉法 ?? ?=+=+),(111n n n n y x hf k k y y 每步计算f 的值一次，其截断误差为O(2 h )。对于预估—校正法。 ()() ???? ?? ? ++==++=+1 212 11 ,,2121k y h x hf k y x hf k k k y y n n n n n n 每步计算f 的值两次，其截断误差为 O(3 h ) 。 下面对预估—校正法进行改进，将该公式写成更一般的形式 ()() bh y ah x hf k y x hf k k R k R y y n n n n n n ++==++=+,,2122111 (3) 其中b a R R ,,,21为待定常数。选择这 些常数的原则是在)(n n x y y =的前提下，使 11)(++-n n y x y ) 的阶尽量 高。为此，作泰勒展开 +++=+++=++=)(] [),(2 12y x y x n n bff af h hf f bk ahf f h bk y ah x hf k 其中 ,,,y x f f f 都是在),(n n y x 处的函数值。将21,k k 代入1+n y 得 +++'++=+++++=+x x n n y x n n ff bR f aR h x y R R h x y ff bR f aR h hf R R y y 22221222 211()()()()()( 与泰勒展开式(2)进行比较，要使得 ()()3 11h o y x y n n =-++ 只要 四个参数满足 ??? ? ? ???? = ==+2 1211 2 221bR aR R R (4) 若1,2 121=== =b a R R , 即得预估—校正公式。 满足(4)式的b a R R ,,,21可以 有各种不同的取法，但不管如何取法，都要计算两次f 的值(即 计算f 在两个不同点的函数值)，截断误差都是()3 h O 。满足条件(4)的一族公式(3)统称为 二 阶龙格—库塔公式。 人们容易想到，如果不增加计算函数值的次数，能否适当地选择这 四个参数，使近似公式的局部截断误差的阶再提高，比如达到 ()4 h O 。为此，把2 k 多展 开 一项，有 ()()() 4 2 22 3 2 222 h O f f b abff f a h bff af h hf k yy xy xx y x ++++ ++= （4）所 以 ()()()()()() 4 2 22222 3 22 2 21122 h O f f R b ff abR f R a h ff bR f aR h x y R R h x y y yy xy x y x n n n ++++ ++'++=+ 而()1+n x y 在n x 的泰勒展开式 ()()()()()() ()()()() 4 2 2 3 2 4 3 2 126 2 6 2 h O ff f f f f ff f h ff f h hf x y h O x y h x y h x y h x y x y y y x yy xy xx y x n n n n n n ++++++ ++ +=+'''+''+'+=+ 由于 () 1+n x y 展开式的3 h 项中 2 y y x ff f f +是无法通过选择参数b a R R ,,,21来消 去的，所以不论四个参 数如何选择，都不能使局部截断误差 ()11++-n n y x y 达到()4 h O 。要想提 高近似公式的阶，只能继续增加计算f 的值的次数。 如 果每步计算三次f 的值，可将公式写 成下列形式： () ()()??? ?? ? ?++=++==+++=+2321313312122 13322111,,,,k b k b y h a x hf k k b y h a x hf k y x hf k k R k R k R y y n n n n n n n n (5) 与二阶龙格—库塔公式的讨论方法类似，要使 ()()4 11h O y x y n n =-++， 只 需8个参数满足?? ? ??? ????? ??== +=++===++61312 113322323222332 232313212 321R b a R a R a R a R a b b a b a R R R 方程组包含6个方程，8个未知量，其解不唯一。满足条件(6)的一族公式(5)统称为三阶龙格—库 塔公式。 一个比较简单的三阶龙格—龙塔公式是 ()()?? ???????+-+=? ?? ?? ++==+++=+2131213 211 2,21,21,6 16461k k y h x hf k k y h x hf k y x hf k R R R y y n n n n n n n n 截断误差为()5 h O 的四阶龙格—库塔公式 是常用的公式，每步都要计算四次 f 的值。它的一般形式是 () () ????? ?? ??++++=+++=++==++=+343242141442 3213133121221443322111,,),(),(k b k b k b y h a x hf k k b k b y h a x hf k k b y h a x h k y x hf k k R k R k R k R y y n n n n n n n n n n (6) (6)式中13个待定常数需满足下列11个方程的方程组 ()()()?? ????? ????? ???? ?? ??? ??= =++=++= ++=++= ++=++++=+===+++241 8 1 1216141 31211443322433422443 323243234222433222433422433224 3433323242 432322244332243 4241432313212 4321R b b a b a b a R a R b a a b a b a R R b a b a b a R R b a R a R a R a R a R a R a R a R a R a b b b a b b a b a R R R R 最常用的四阶龙格—库 塔公式是标准四阶 龙格—库塔公式： ()() () ???? ?? ????? ++=?? ? ??++=??? ??++==++++=+342312 1 43211 ,21,2121,21,226 1k y h x hf k k y h x hf k k y h x hf k y x hf k k k k k y y n n n n n n n n n n 和吉尔公式 ()??? ??? ???? ?? ?? ??? ????? ?? ++??? ??-++=??? ? ???? ? ??-+??? ??+-++=??? ??++==+??? ??++??? ??-++=+3242231214 32 1121121,2112121,2121,21,61211312113161k k y h x hf k k k y h x hf k k y h x hf k y x hf k k k k k y y n n n n n n n n n n 标准四阶龙格—库塔公式手算时采用表8-2 所示的表格计算。 例：用标准四阶龙格—库塔方法求解初值问题 ??? ? ? =≤≤=-=2 .0,6.00,1)0(2h x y y x y dx dy 解：计算过程和结果 如表7-3所示 因此有 表7-3 ()()()()483281 .16.0,341667.14.0183229 .12.0,10≈≈≈=y y y y 对该例，用几种不同的一步法计算的结果如表7-4。 由表7-4可见，虽然四阶龙格—库塔方法每步要计算四次f 的值，但以h =0.2 为 步长的计算结果就有5位有效数字，而欧拉法与预估—校正法以h =0.1 为步长的计算结果 才具有2位与3位有效数字。如果步长h 也取0.2，则结果的精度会更低。 表7-4 求解初值问题的一步法在计算时只用到前面一步的结果，所以要提高精度时，需要增加中间函数值的计算，这就加大了计算量。如果在计算y n+1时，不仅用到x n 上的近似值y n ，还用到前面若干节点x n-1 ,x n-2 ，…上的近似值y n-1 ,y n-2 ,…，称这种方法为多步法。亚当斯(Adams)方法是多步法中的一种。 我们知道，求解初值问题 0)() ,(y x y y x f dx dy == (1) 等价于求解积分方程 ? ++ =+1 ))(,()()(1n n x x n n dx x y x f x y x y (2) 选用不同的数值方法计算(8.3.2)式右端的积分项，就会导出不同的计算公式。例如，用矩形法计算积分项 )) (,())(,(1 n n x x x y x hf dx x y x f n n ≈? + 代入(2)式得 )) (,()()(1n n n n x y x hf x y x y +=+ 离散化即得欧拉公式，其截断误差为O(h 2)。 为了提高精度，改用梯形法计算积分项 ))] (,())(,([2 ))(,(111 +++≈ ? +n n n n x x x y x f x y x f h dx x y x f n n 代入(2)式得 ))] (,())(,([2 )()(111+++++ =n n n n n n x y x f x y x f h x y x y 离散化得到改进的欧拉公式，其截断误差是O(h 3)。 为此，我们想到，基于插值原理可以建立一系列的数值积分方法，运用这些方法可以导出求初值问题的一系列计算公式。一般地，若用插值多项式ψk (x)，代替 f(x,y(x)),用?+1)(n n x x k dx x ?作为 ? +1 ))(,(n n x x dx x y x f 的近似值，离散化后 得公式 ? ++ =+1 )(1n n x x k n n dx x y y ? (3) 7.3.1亚当斯显式公式 设初值问题(1)的解y(x)在x n ,x n-1 ,…,x n-k 上各点的近似值y n ,y n-1 ,…,y n-k 都已计算出来，用点(x n ,f n ),(x n-1 ,f n-1),…，(x n-k ,f n-k )作后插公式 ∑ =-++?= -++? + ++?+ ??==k j n j K n n n k j t t t j f k t t t K t t f t f f x 0 2 ) 1()1(! ) 1()1(! )1(! 2)( ? 其中x=x n +th ，将ψk (x)代入(3)式，得亚当斯显式公式 ∑=+?+=k j n j j n n f a h y y 0 1 其中 ? -++= 1 )1()1(! 1dt j t t t j a j 它的前几个值见表7-5． 表7-5 当 k=0 时，(8.3.4)式为欧拉公式。 当 k=1 时，(8.3.4)式为 ) 21(1n n n n f f h y y ?+ +=+ 当 k=3 时，(4)式为 ) 8 312 52 1(3 2 1n n n n n n f f f f h y y ?+ ?+ ?+ +=+ (5) 利用差分与函数值之间的关系，(5)可写成下列形式 ) 9375955(24 3211---+-+-+ =n n n n n n f f f f h y y (6) 假设(6) 式右端的值精确，即y n-i =y(x n-i )，i=0,1,2,3. 那么 f n-i =f(x n-i ,y n-i )=f(x n-i ,y(x n-i ))=y ′(x n-i ) 代入(8.3.6)式 )) (9)(37)(59)(55(24 3211---+'-'+'-'+ =n n n n n n x y x y x y x y h y y 将上式右端各项在点x n 处展开，得 +- + '''+ ''+ '+=+)(144 49)(24 )(6 )(2 )()() 5(5) 4(4 3 2 1n n n n n n n x y h x y h x y h x y h x y h x y y 另一方 面，对于准确解y(x n+1)，有 +- + '''+ ''+ '+=+)(120 )(24 )(6 )(2 )()()() 5(5 ) 4(4 3 2 1n n n n n n n x y h x y h x y h x y h x y h x y x y 于是 ) ()(720 251)(5 ) 5(5 11h O x y h y x y n n n =+= -++ 所以公式(6)的局部截断误差为O(h 5)，它是四阶方法，称(6)式为四阶亚当斯显式公式。 7.3.2亚当斯隐式公式 亚当斯显式方法选取x n ,x n-1 ,…,x n-k 作为插值点，这时，用插值函数ψk (x)在求积区间［x n ,x n+1 ］上逼近函数 f(x,y(x)) 实际上是一个外推过程，效果不太理想。为了提高精度，我们变外推为内插，改用 x n+1 ,x n ,…,x n-k+1为插值节点建立后插公式ψk (x),然后重复前一段的推导过程，得到亚当斯隐式公式 ∑=++?+=k j n j j n n f b h y y 0 1 1 其中 ?--++= 1 )1()1(! 1 dt j t t t j b j 它的前几个值见表7-6。 表7-6 当 k=0 时，(7)式为 11+++=n n n hf y y 当 k=1 时，(7)式为梯形公式 当 k=3 时，(7)式为 ) 24 112 12 1(13 12 111+++++?- ?- ?- +=n n n n n n f f f f h y y 即 ) 15199(24 2111++-++ =-+++n n n n n n f f f f h y y (8) 类似于亚当斯显式公式，可以导出隐式公式(8.3.8)的局部截断误差为 ) ()(720 19)(5 ) 5(5 11h O x y h y x y n n n =+- =-++ 所以隐式公式(8)是四阶方法。 7.3.3亚当斯预估一校正公式 把亚当斯显式公式与隐式公式联立使用，构造亚当斯预估一校正公式。以四阶亚当斯方法为例有下列公式： ) f 9f 37f 59f 55(24 h y y 3n 2n 1n n n ) 0(1n ---+-+-+ = ) f f 5f 19y ,x (f 9(24 h y y 2n 1n n ) k (1n 1n n ) 1k (1 n --+++++-++= (9) 迭代过程进行到连续 两次迭代结果之差的绝对值小于给定的精度ε为止。 这时取) 1k (1n 1n y y +++=，然后转入下一步计算。 只迭代1次的公式 ) f 9f 37f 59f 55(24h y y 3n 2n 1n n n ) 0(1n ---+-+-+ = ) f f 5f 19y ,x (f 9(24 h y y 2n 1n n ) k (1n 1n n ) 1k (1 n --+++++-++ = (.10) 称为亚当斯预估一校正公式。第1式称为预估公式，第2式称为校正公式。其局部截断误差为 ) h (O )x (y h 720 19y )x (y 5 n ) 5(51n 1n =- ≈-++ 下面我们具体讨论误差估计问题。设 由(10)式求得准确解 y(x n+1)的近似解为y n+1，由于 ) x (y h 720251y )x (y n ) 5(5) 0(1n 1n ≈ -++ ) x (y h 720 19y )x (y n ) 5(51n 1n - ≈-++ 上面二式相减得 ) x (y h 720 270y )x (y n ) 5(5) 0(1n 1n ≈ -++ 所以 27019y y y )x (y ) 0(1 n 1n 1n 1n - ≈--++++ 于是有 ) y y (270 19y )x (y ) 0(1n 1n 1n 1n ++++-- ≈- 若 <--++) 0(1n 1n y y (270 19ε 则可持续求解y n+2,…，否则,应缩小步长h 重新计算。 与同阶的龙格—库塔方法相比较，亚当斯方法计算量小，公式简单，程序易于实现。但是由于亚当斯方法在计算y n+1时，不仅用到前面一步的信息，而且还要用到更前面几步的信息，因此它不能自动开始。开始的前几个值必须借助于一步法获得。 例：分别用四阶亚当斯显式公式与预估一校正公式求解初值问题 y x 2y dy dx - = 1)0(y = 1.0h ,1x 0=≤≤ 解：用第二节例1，按标准四阶龙格—库塔方法求得的结果y 1,y 2,y 3作为开始值，然后用显式公式(6)与预估一校正方法(10) 进行计算，计算结果如表7-7所示。 表7-7 近似解与准确解进行比较，显式方法的结果有4位有效数字，而预估一校正法的结果则有7位有效数字。 第四节 线性多步法 线性多步法的一般公式为 1+n y =n y 0α+ 1 1-n y α+…+k n k y -α+h(k n k n n f f f -+-+++βββ...011) (1) 当β-1 =0时，公式(1)是显式的；当β-1 ≠0时 ，公式(1)是隐式的。“线性”是指 1+n y 是n y ,1-n y , … , k n y -, 1+n f ,n f ,…,k n f -的线性组 合，“多步”是指计算y n+1时 ，不仅用到前一步的结果y n ,而且用到更前面几步的结果。 所有形如(8.4.1)式的公式都 可以利用泰勒展开的方法构造出来。作为例子，我们推导形如 1+n y = n y 0α+ 11-n y α+ 2 2-n y α+h( 3 322110---+++n n n n f f f f ββββ) (2) 的数值计算公式。 假设 i y =y(i x ), i=n, n-1, n-2 i f =f(i x ,i y ) i=n,n-1,n-2,n-3 将上述各项在 x n 处泰勒展开 Y n-i =y(x n -ih)=y(x n )-ih ) (' n x y + )(! 2)(' '2 n x y ih - )(! 3)(' ''3 n x y ih +…, i=1,2 i n f -= )(' i n x y -= ) (' ih x y n - =) (' n x y - ) (n x y ih '+ )(! 2)(2 n x y ih '''- )(! 3)() 4(3 n x y ih +…, i=1,2,3 代入(2)式 得 1 +n y =(210ααα++)y(n x )+(3210202ββββαα++++--))(n x y h '+ (321216424βββαα---+))(! 22 n x y h ''+ ( 3 2121271238βββαα+++--) )(! 33 n x y h '''+ ( 3 212110832416βββαα---+) )(! 4) 4(4 n x y h + (3) 将 y(1+n x )在n x 处泰勒展开 ...)(! 3)(! 3)(! 2)()()(3 3 2 1+'''+ '''+ ''+ '+=+n n n n n n x y h x y h x y h x y h x y x y (4) 比较(3)式与(4)式，让h 的同次幂的系数相等，得 1210=++ααα 12321021=++++--ββββαα 1642432121=---+βββαα (5) 127123832121=+++--βββαα 11083241632121=---+βββαα …… 如果(5)的前7个方程求出未知数3210210,,,,,,ββββααα,则其局部截断误差为 )(7 h O .当然，也可以从前k 个方程解出未知数(k ≤7)，其局部截断误差为 )(k h O .比如从前5 个方程解出 7个未知数，因此，有2个自由未知数，若令 021==αα,则从(5)的前5 个方程解得 24 9,24 37,2459,2455,132100- == - == =ββββα 对应的公式是 ) 9375955(24 3211---+-+-+ =n n n n n n f f f f h y y 它正好是四阶阿达姆斯显式公式。 类似地， 若设计算公式为 1 +n y = n y 0α+11-n y α+22-n y α+h(2211011--+-+++n n n n f f f f ββββ) 为使这类公式为四阶方法，系数应满足 1210=++ααα 12210121=++++---ββββαα 1422421121=--++-βββαα (7) 11283821121=+++---βββαα 132441621121=--++-βββαα 特别地，取021==αα，则从(7)式解得 24 1,245,2419,24 9,121010= - =- == =-ββββα 对应的公式是 ) 5199(24 2111--+++-++=n n n n n n f f f f h y y 它正好是四 阶亚当斯隐式公式。 第五节 方程组与高阶方程的数值解法 前面几节介绍了一阶微分方程初值问题的各种数值解法，这些解法同样适用 于一阶微分方程组与高阶方程。为了避免书写的复杂，下面仅就两个未知函数的方程组和二阶方程为例说明各种方法的计算公式。 7.5.1一阶方程组 设有一阶微分方程组初值问题 第8章 常微分方程边值问题的数值解法 引 言 第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。 只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为 则边值问题(8.1.1)有唯一解。 推论 若线性边值问题 ()()()()()(),, (),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤?? ==? (8.1.2) 满足 (1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。 求边值问题的近似解，有三类基本方法： (1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解； (2) 有限元法(finite element method)； (3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。 差分法 8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=< 第 一 章 一阶微分方程的解法的小结 ⑴、可分离变量的方程： ①、形如 )()(y g x f dx dy = 当0)(≠y g 时，得到 dx x f y g dy )() (=，两边积分即可得到结果； 当0)(0=ηg 时，则0)(η=x y 也是方程的解。 例1.1、 xy dx dy = 解：当0≠y 时，有 xdx y dy =，两边积分得到)(2ln 2为常数C C x y += 所以)(112 12 C x e C C e C y ±==为非零常数且 0=y 显然是原方程的解； 综上所述，原方程的解为)(12 12 为常数C e C y x = ②、形如0)()()()(=+dy y Q x P dx y N x M 当0)()(≠y N x P 时，可有 dy y N y Q dx x P x M ) () ()()(=，两边积分可得结果； 当0)(0=y N 时，0y y =为原方程的解，当0(0=） x P 时，0x x =为原方程的解。 例1.2、0)1()1(2 2 =-+-dy x y dx y x 解：当0)1)(1(2 2 ≠--y x 时，有 dx x x dy y y 1 122-=-两边积分得到 )0(ln 1ln 1ln 22≠=-+-C C y x ，所以有)0()1)(1(22≠=--C C y x ； 当0)1)(1(2 2 =--y x 时，也是原方程的解； 综上所述，原方程的解为)()1)(1(2 2 为常数C C y x =--。 ⑵可化为变量可分离方程的方程： ①、形如 )(x y g dx dy = 解法：令x y u =，则udx xdu dy +=，代入得到)(u g u dx du x =+为变量可分离方程，得到 常微分方程解题方法总结 来源：文都教育 复习过半， 课本上的知识点相信大部分考生已经学习过一遍 . 接下来， 如何将零散的知 识点有机地结合起来， 而不容易遗忘是大多数考生面临的问题 . 为了加强记忆， 使知识自成 体系，建议将知识点进行分类系统总结 . 著名数学家华罗庚的读书方法值得借鉴， 他强调读 书要“由薄到厚、由厚到薄”，对同学们的复习尤为重要 . 以常微分方程为例， 本部分内容涉及可分离变量、 一阶齐次、 一阶非齐次、 全微分方程、 高阶线性微分方程等内容， 在看完这部分内容会发现要掌握的解题方法太多， 遇到具体的题 目不知该如何下手， 这种情况往往是因为没有很好地总结和归纳解题方法 . 下面以表格的形 式将常微分方程中的解题方法加以总结，一目了然，便于记忆和查询 . 常微分方程 通解公式或解法 ( 名称、形式 ) 当 g( y) 0 时，得到 dy f (x)dx ， g( y) 可分离变量的方程 dy f ( x) g( y) 两边积分即可得到结果； dx 当 g( 0 ) 0 时，则 y( x) 0 也是方程的 解 . 解法：令 u y xdu udx ，代入 ，则 dy 齐次微分方程 dy g( y ) x dx x u g (u) 化为可分离变量方程 得到 x du dx 一 阶 线 性 微 分 方 程 P ( x)dx P ( x) dx dy Q(x) y ( e Q( x)dx C )e P( x) y dx 伯努利方程 解法：令 u y1 n，有 du (1 n) y n dy ， dy P( x) y Q( x) y n（n≠0,1）代入得到du (1 n) P(x)u (1 n)Q(x) dx dx 求解特征方程：2 pq 三种情况： 二阶常系数齐次线性微分方程 y p x y q x y0 二阶常系数非齐次线性微分方程 y p x y q x y f ( x) (1)两个不等实根：1, 2 通解： y c1 e 1x c2 e 2x (2) 两个相等实根：1 2 通解： y c1 c2 x e x (3) 一对共轭复根：i , 通解： y e x c1 cos x c2 sin x 通解为 y p x y q x y 0 的通解与 y p x y q x y f ( x) 的特解之和. 常见的 f (x) 有两种情况： x （ 1）f ( x)e P m ( x) 若不是特征方程的根，令特解 y Q m ( x)e x；若是特征方程的单根，令特 解 y xQ m ( x)e x；若是特征方程的重根， 令特解 y*x2Q m (x)e x； （2）f (x) e x[ P m ( x) cos x p n ( x)sin x] 淮北师范大学 2013届学士学位论文 常微分方程数值解法的误差分析 学院、专业数学科学学院数学与应用数学 研究方向计算数学 学生姓名李娜 学号 20091101070 指导教师姓名陈昊 指导教师职称讲师 年月日 常微分方程数值解法的误差分析 李娜 （淮北师范大学数学科学学院，淮北，235000） 摘要 自然界与工程技术中的很多现象，往往归结为常微分方程定解问题。许多偏微分方程问题也可以化为常微分方程问题来近似求解。因此，研究常微分方程的数值解法是有实际应用意义的。数值解法是一种离散化的数学方法，可以求出函数的精确解在自变量一系列离散点处的近似值。随着计算机计算能力的增强以及数值计算方法的发展，常微分方程的数值求解方法越来越多，比较成熟的有Euler 法、后退Euler法、梯形方法、Runge—Kutta方法、投影法和多步法，等等．本文将对这些解的误差进行分析,以求能够得到求解常微分数值解的精度更好的方法。 关键词：常微分方程, 数值解法, 单步法, 线性多步法, 局部截断误差 Error Analysis of Numerical Method for Solving the Ordinary Differential Equation Li Na (School of Mathematical Science, Huaibei Normal University, Huaibei, 235000) Abstract In nature and engineering have many phenomena , definite solution of the problem often boils down to ordinary differential equations. So study the numerical solution of ordinary differential equations is practical significance. The numerical method is a discrete mathematical methods, and exact solution of the function can be obtained in the approximation of a series of discrete points of the argument.With the enhanced computing power and the development of numerical methods,ordinary differential equations have more and more numerical solution,there are some mature methods. Such as Euler method, backward Euler method, trapezoidal method, Runge-Kutta method, projection method and multi-step method and so on.Therefore, numerical solution of differential equation is of great practical significance. Through this paper, error of these solutions will be analyzed in order to get a the accuracy better way to solve the numerical solution of ordinary differential. Keywords:Ordinary differential equations, numerical solution methods, s ingle ste p methods, l inear multi-step methods, local truncation error i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。 常微分方程初值问题数值解法 朱欲辉 （浙江海洋学院数理信息学院, 浙江舟山316004） [摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实 际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正 公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点. [关键词]：常微分方程;初值问题; 数值方法; 收敛性; 稳定性; 误差估计 Numerical Method for Initial-Value Problems Zhu Yuhui (School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004) [Abstract]:In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be e xpressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. [Keywords]:Ordinary differential equation; Initial-value problem; Numerical method; Convergence; Stability;Error estimate 各类微分方程的解法 1.可分离变量的微分方程解法 一般形式:g(y)dy=f(x)dx 直接解得∫g(y)dy=∫f(x)dx 设g(y)及f(x)的原函数依次为G(y)及F(x),则G(y)=F(x)+C为微分方程的隐式通解 2.齐次方程解法 一般形式:dy/dx=φ(y/x) 令u=y/x则y=xu,dy/dx=u+xdu/dx,所以u+xdu/dx=φ(u),即du/［φ(u)-u］=dx/x 两端积分,得∫du/［φ(u)-u］=∫dx/x 最后用y/x代替u,便得所给齐次方程的通解 3.一阶线性微分方程解法 一般形式:dy/dx+P(x)y=Q(x) 先令Q(x)=0则dy/dx+P(x)y=0解得y=Ce－ ∫P(x)dx,再令y=u e－∫P(x)dx代入原方程解得u=∫Q(x) e∫P(x)dx dx+C,所以y=e－∫P(x)dx［∫Q(x)e∫P(x)dx dx+C］ 即y=Ce－∫P(x)dx +e－ ∫P(x)dx∫Q(x)e∫P(x)dx dx为一阶线性微分方程的通解 4.可降阶的高阶微分方程解法 ①y(n)=f(x)型的微分方程 y(n)=f(x) y(n-1)= ∫f(x)dx+C1 y(n-2)= ∫［∫f(x)dx+C1］dx+C2 依次类推,接连积分n次,便得方程y(n)=f(x)的含有n个任意常数的通解②y”=f(x,y’) 型的微分方程 令y’=p则y”=p’,所以p’=f(x,p),再求解得p=φ(x,C1) 即dy/dx=φ(x,C1),所以y=∫φ(x,C1)dx+C2 ③y”=f(y,y’) 型的微分方程 令y ’=p 则y ”=pdp/dy,所以pdp/dy=f(y,p),再求解得p=φ(y,C 1) 即dy/dx=φ(y,C 1),即dy/φ(y,C 1)=dx,所以∫dy/φ(y,C 1)=x+C 2 5.二阶常系数齐次线性微分方程解法 一般形式:y ”+py ’+qy=0,特征方程r 2+pr+q=0 6.二阶常系数非齐次线性微分方程解法 一般形式: y ”+py ’+qy=f(x) 先求y ”+py ’+qy=0的通解y 0(x),再求y ”+py ’+qy=f(x)的一个特解y*(x) 则y(x)=y 0(x)+y*(x)即为微分方程y ”+py ’+qy=f(x)的通解 求y ”+py ’+qy=f(x)特解的方法: ① f(x)=P m (x)e λx 型 令y*=x k Q m (x)e λx ［k 按λ不是特征方程的根,是特征方程的单根或特征方程的重根依次取0,1或2］再代入原方程,确定Q m (x)的m+1个系数 ② f(x)=e λx ［P ｌ(x)cos ωx+P n (x)sin ωx ］型 令y*=x k e λx ［Q m (x)cos ωx+R m (x)sin ωx ］［m=max ﹛ｌ,n ﹜,k 按λ+i ω不是特征方程的根或是特征方程的单根依次取0或1］再代入原方程,分别确定Q m (x)和R m (x)的m+1个系数 第八章 常微分方程数值解法 考核知识点： 欧拉法，改进欧拉法，龙格-库塔法，单步法的收敛性与稳定性。 考核要求： 1． 解欧拉法，改进欧拉法的基本思想；熟练掌握用欧拉法，改进欧拉法、求微 分方程近似解的方法。 2． 了解龙格-库塔法的基本思想；掌握用龙格-库塔法求微分方程近似解的方 法。 3． 了解单步法的收敛性、稳定性与绝对稳定性。 例1 用欧拉法，预估——校正法求一阶微分方程初值问题 ? ??=-='1)0(y y x y ，在0=x （0，1）0.2近似解 解 （1）用1.0=h 欧拉法计算公式 n n n n n n x y y x y y 1.09.0)(1.01+=-+=+，1.0=n 计算得 9.01=y 82.01.01.09.09.02=?+?=y （2）用预估——校正法计算公式 1,0)(05.01.09.0)0(111)0(1=???-+-+=+=++++n y x y x y y x y y n n n n n n n n n 计算得 91.01=y ，83805.02=y 例2 已知一阶初值问题 ???=-='1 )0(5y y y 求使欧拉法绝对稳定的步长n 值。 解 由欧拉法公式 n n n n y h y h y y )51(51-=-=+ n n y h y ~)51(~1-=+ 相减得01)51()51(e h e h e n n n -==-=-Λ 当 151≤-h 时，4.00≤ 第八章 常微分方程的数值解法 一．内容要点 考虑一阶常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节 点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。 用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。 (一)常微分方程处置问题解得存在唯一性定理 对于常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 如果： (1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。 (2) ),(y x f 对于y 满足利普希茨条件，即 2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程?????==0 0)() ,(y x y y x f dx dy 的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。 定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。 收敛性定理：若一步方法满足： (1)是p 解的. (2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件. (3) 初始值y 0是精确的。则),()()(p h O x y kh y =-kh =x -x 0,也就是有 0x y y lim k x x kh 0h 0 =--=→)( (一)、主要算法 1．局部截断误差 局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~ +k y 的误差y (x k+1)- 1~ +k y 称为局部截断误差。 注意：y k+1和1~ +k y 的区别。因而局部截断误差与误差e k +1=y (x k +1) -y k +1不同。 如果局部截断误差是O (h p+1),我们就说该数值方法具有p 阶精度。 常微分方程数值解法 【作用】微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步： 1. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。 2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。 3. 运用这些规律列出方程和定解条件。基本模型 1. 发射卫星为什么用三级火箭 2. 人口模型 3. 战争模型 4. 放射性废料的处理通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来” 的于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 1. 改进Euler 法： 2. 龙格—库塔( Runge—Kutta )方法： 【源程序】 1. 改进Euler 法： function [x,y]=eulerpro(fun,x0,x1,y0,n);%fun 为函数，(xO, x1)为x 区间，yO 为初始值，n 为子 区间个数 if nargin<5,n=5O;end h=(x1-xO)/n; x(1)=xO;y(1)=yO; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(fun,x(i),y(i)); y2=y(i)+h*feval(fun,x(i+1),y1); y(i+1)=(y1+y2)/2; end 调用command 窗口 f=i nlin e('-2*y+2*x A2+2*x') [x,y]=eulerpro(f,O,,1,1O) 2 x +2x , (0 < x < , y(0) = 1 求解函数y'=-2y+2 2. 龙格—库塔( Runge—Kutta )方法： [t,y]=solver('F',tspan ，y0) 这里solver为ode45, ode23, ode113，输入参数F是用M文件定义的微分方程y'= f (x, y)右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。 注：ode45和ode23变步长的，采用Runge-Kutta算法。 ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长(变步长)的常微分方程数值解法，其整体截断误差为(△ 口人5解 决的是Nonstiff(非刚性)常微分方程。 实验4 常微分方程的数值解 【实验目的】 1．掌握用MATLAB软件求微分方程初值问题数值解的方法； 2．通过实例用微分方程模型解决简化的实际问题； 3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。 在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式： a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8 dh/dt=v 又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下： function [ dx ] = rocket( t,x ) a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8; dx=[x(2);a]; end ts=0:1:60; x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a]; 数据如下： t h v a 0 0 0 13.06 1.00 6.57 13.19 13.30 2.00 26.44 26.58 1 3.45 3.00 59.76 40.06 13.50 4.00 106.57 53.54 13.43 5.00 16 6.79 66.89 13.26 6.00 240.27 80.02 12.99 7.00 326.72 92.83 12.61 8.00 425.79 105.22 12.15 9.00 536.99 117.11 11.62 10.00 659.80 128.43 11.02 11.00 793.63 139.14 10.38 12.00 937.85 149.18 9.71 13.00 1091.79 158.55 9.02 14.00 1254.71 167.23 8.33 15.00 1425.93 175.22 7.65 16.00 1604.83 182.55 6.99 17.00 1790.78 189.22 6.36 18.00 1983.13 195.27 5.76 19.00 2181.24 200.75 5.21 20.00 2384.47 205.70 4.69 21.00 2592.36 210.18 4.22 22.00 2804.52 214.19 3.79 23.00 3020.56 217.79 3.41 24.00 3240.08 221.01 3.07 25.00 3462.65 223.92 2.77 26.00 3687.88 226.56 2.50 27.00 3915.58 228.97 2.27 -293- 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如 22x y dx dy +=。于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 §1 常微分方程的离散化 下面主要讨论一阶常微分方程的初值问题，其一般形式是 ?????=≤≤=0 )() ,(y a y b x a y x f dx dy （1） 在下面的讨论中，我们总假定函数),(y x f 连续，且关于y 满足李普希兹(Lipschitz)条 件，即存在常数L ，使得 |||),(),(|y y L y x f y x f ?≤? 这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。 所谓数值解法，就是求问题（1）的解)(x y 在若干点 b x x x x a N =<<<<=L 210 处的近似值),,2,1(N n y n L =的方法，),,2,1(N n y n L =称为问题（1）的数值解， n n n x x h ?=+1称为由n x 到1+n x 的步长。今后如无特别说明，我们总取步长为常量h 。 建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i ）用差商近似导数 若用向前差商h x y x y n n ) ()(1?+代替)('n x y 代入（1）中的微分方程，则得 )1,,1,0())(,() ()(1?=≈?+N n x y x f h x y x y n n n n L 化简得 ))(,()()(1n n n n x y x hf x y x y +≈+ 如果用)(n x y 的近似值n y 代入上式右端，所得结果作为)(1+n x y 的近似值，记为1+n y ， 则有 )1,,1,0() ,(1?=+=+N n y x hf y y n n n n L （2） 这样，问题（1）的近似解可通过求解下述问题 ?? ?=?=+=+) () 1,,1,0(),(01a y y N n y x hf y y n n n n L （3） 得到，按式（3）由初值0y 可逐次算出N y y y ,,,21L 。式（3）是个离散化的问题，称为差分方程初值问题。 i.常微分方程初值问题数值解法 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法--差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<<<<=L （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-L 方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t L 上的差分解1,,N u u L 。 二阶微分方程： 时为非齐次 时为齐次，0)(0)()()()(22≠≡=++x f x f x f y x Q dx dy x P dx y d 二阶常系数齐次线性微分方程及其解法： 2 122,)(2,,(*)0)(1,0(*)r r y y y r r q pr r q p qy y p y 式的两个根、求出的系数； 式中的系数及常数项恰好是，，其中、写出特征方程：求解步骤： 为常数； ，其中?'''=++?=+'+''式的通解：出的不同情况，按下表写、根据(*),321r r 二阶常系数非齐次线性微分方程 型 为常数； 型，为常数 ，]sin )(cos )([)()()(,)(x x P x x P e x f x P e x f q p x f qy y p y n l x m x ωωλλλ+===+'+'' 二阶常系数非齐次线性微分方程的一般形式是 ''+'+=y py qy f x () （1） 其中p q ,是常数。 方程（1）的通解为对应的齐次方程 0=+'+''qy y p y （2） 的通解Y 和方程（1）的一个特解*y 之和。即 *y Y y +=.我们已解决了求二阶常系数齐 次线性方程通解的问题，所以，我们只需讨论求二阶常系数非齐次线性微分方程的特解* y 的方法。 下面我们只介绍当方程（1）中的)(x f 为如下两种常见形式时求其特解*y 的方法。 一、 f x e P x x m ()()=?λ型 由于方程（1）右端函数f x ()是指数函数e x λ?与m 次多项式P x m ()的乘积，而指数 函数与多项式的乘积的导数仍是这类函数，因此，我们推测： 方程（1）的特解应为 y e Q x x *?=λ()( Q x ()是某个次数待定的多项式 ) y e Q x e Q x x x *??'=+'λλλ()() y e Q x Q x Q x x *?"=?+'+''λλλ[()()()]22 代入方程（1），得 e Q x p Q x p q Q x e P x x x m λλλλλ???''++'+++≡?[()()()()()]()22 消去e x λ?，得 ''++'+++≡Q x p Q x p q Q x P x m ()()()()()()22λλλ （3） 讨论 01、如果λ不是特征方程 r pr q 20++=的根。 即 02≠++q p λλ 由于P x m ()是一个m 次的多项式，欲使（3）的两端恒等，那未Q x ()必为一个m 次多项式，设为 Q x b x b x b x b m m m m m ()=++++--0111Λ 将之代入（3），比较恒等式两端x 的同次幂的系数，就得到以b b b b m m 01 1,,,,Λ-为未知数的m +1个线性方程的联立方程组，解此方程组可得到这m +1个待定的系数，并得到特解 y e Q x x m *?=λ() 02、如果λ是特征方程 r pr q 20++=的单根。 即 λλ20++=p q ，但 20λ+≠p 欲使（3）式的两端恒等，那么'Q x ()必是一个m 次多项式。 因此，可令 Q x x Q x m ()()=? 并且用同样的方法来确定)(x Q 的系数b b b b m m 0 11,,,,Λ-。 03、如果λ是特征方程 r pr q 20++=的二重根。 即 λλ20++=p q ，且 20λ+=p 。 欲使（3）式的两端恒等，那么''Q x ()必是一个m 次多项式 因此， 可令 Q x x Q x m ()()=?2 并且用同样的方法来确定)(x Q 的系数b b b b m m 011,,,,Λ-。 目录 1 引言 (1) 2 二阶常系数常微分方程的几种解法 (1) 2.1 特征方程法 (1) 2.1.1 特征根是两个实根的情形 (2) 2.1.2 特征根有重根的情形 (2) 2.2 常数变异法 (4) 2.3 拉普拉斯变化法 (5) 3 常微分方程的简单应用 (6) 3.1 特征方程法 (7) 3.2 常数变异法 (9) 3.3 拉普拉斯变化法 (10) 4 总结及意义 (11) 参考文献 (12) 二阶常微分方程的解法及其应用 摘要：本文通过对特征方程法、常数变易法、拉普拉斯变换法这三种二阶常系数常微分方程解法进行介绍，特别是其中的特征方程法分为特征根是两个实根的情形和特征根有重根的情形这两种情况，分别使用特征值法、常数变异法以及拉普拉斯变换法来求动力学方程，现今对于二阶常微分方程解法的研究已经取得了不少成就，尤其在二阶常系数线性微分方程的求解问题方面卓有成效。应用常微分方程理论已经取得了很大的成就，但是，它的现有理论也还远远不能满足需要，还有待于进一步的发展，使这门学科的理论更加完善。 关键词：二阶常微分方程;特征分析法；常数变异法；拉普拉斯变换 METHODS FOR TWO ORDER ORDINARY DIFFERENTIAL EQUATION AND ITS APPLICATION Abstract:This paper introduces the solution of the characteristic equation method, the method of variation of parameters, the Laplasse transform method the three kind of two order ordinary differential equations with constant coefficients, especially the characteristic equation method which is characteristic of the root is the two of two real roots and characteristics of root root, branch and don't use eigenvalue method, method of variation of constants and Laplasse transform method to obtain the dynamic equation, the current studies on solution of ordinary differential equations of order two has made many achievements, especially in the aspect of solving the problem of two order linear differential equation with constant coefficients very fruitful. Application of the theory of ordinary differential equations has made great achievements, however, the existing theory it is still far from meeting the need, needs further development, to make the discipline theory more perfect. Keywords:second ord er ordinary differential equation; Characteristic analysis; constant variation method; Laplasse transform 1 引言 数学发展的历史告诉我们，300年来数学分析是数学的首要分支，而微分方程常微分方程边值问题的数值解法
一阶常微分方程解法总结
常微分方程解题方法总结.doc
常微分方程数值解法的误差分析教材
常微分方程数值解法
常微分方程初值问题数值解法
各类微分方程的解法大全
常微分方程数值解法
常微分方程数值解法
常微分方程数值解法
常微分方程的数值解
15第十五章 常微分方程的解法
(整理)常微分方程数值解法
(完整版)专题一(二阶常微分方程解法)
二阶常微分方程的解法及其应用
相关主题
文本预览