高阶常微分方程的数值求解'谷照升（长春工程学院理学院，长春，130012）摘要对经典初始条件的高阶常微分方程，给出其数值求解方法。该方法比Runge-Kutta法具冇更好的适应性、易用性、计算速度和可控制的更高精度。关键词常微分方程；数值解：算法中图分类号：0241.81文献标识码:A1引言求解复杂的1阶常微分方程，通常只能采用数值解法。数值解法一般又以Runge-Kutta法为主。对高阶常微分方程，则通常是将其转化为1阶常微分方程组，再用Runge-Kutta法求解⑴％但这种通用性方法，在精度、软件计算的适应度方面，常常不够理想，英至得不到结果。针对不同的广普性方程类型，可以建立更具针对性的计算方法。例如，针对三类边值条件和特定形式的方程，己经有有相应的差分法和有限元法。每一种更具针对性的方法，都有其更髙的楮度和更为健壮的算法，当然也存在其必然的局限性。木文对如下形式的2阶常微分方程^-j-+tz（x）-y-+Z?（x）y=/（X）dx^ax<y（Q）=儿皿S,b］（1）y'（Q）=冗给出了完整、方便、高效的单步法数值求解算法，并将其推广到任意高阶问题的求解。2算法思路：将［a,b］分割为a=xQ<Xi<...<xn=b,步长Ax,=%.-x,.,o从匸1开始，利用数值积分公式,通过［冷旬上的数值积分，先求f（£）,再求出y'a）,最后求yg）。依次取匸2,3,…川得到各点y"（兀J、«/（兀）、y（Q的近似值。根据数值积分的算法不同，主要有两种不同的计算公式。2・1、矩形数值积分算法矩形算法形式简单，粕度偏低。基于不同的积分形式，乂可以得到3种不同的计算公式，其精基金项目：吉林省自然科学基金项H（201215115）度无显著差別。此处仅给出与2.2梯形算法最为接近的一种。利用/（%,）一：/（£•_］）=fy\x）dx=y\xi）^xiJxi-\丿（兀）-y（s）=£y\x）dx«y\xi}\xi在x{点得到y（x,）«y\Xi）Ax；+ya-）⑵yU,）=/（乞）心「+^（^-i）=fa）山「+y（兀-）心,+刃"“）（3）这里i=1,2,...,no将（2）、（3）代入（1）,在兀点得到/（xJ+^Jt/Cx^Ax.+y（G）］+〃（兀J［Axi）Ax：+/（%,_））Ax,.+y（G）］=fg从而/（v）=/a）一aa）y（®j-"兀）2（心|）心,+刃兀―卄（4）z1+a（xi）\xi+b（xj$算法：Step1:z=l;将⑴中y（兀o）=儿、/（xo）=ya代入⑷，输出y"3）;将y"（兀i）代入⑵，输岀y'a）；収y（K）=：/（兀0）心1+y（xo）»输出y（x\）°Step2:whilei<n{z=z+l;将y（兀一i）、y（xi~\）代入（4）,输出y\Xi）；将y"a）代入⑵，输出/（x,）；取y（无）=）/（兀_】）心,+y（£_］）,输出y（兀）；}2.2、梯形数值积分算法根据（1）式，首先算出(5)y\xQ）=/（x0）一«（x0）/（x0）一b（xQ）y（xQ）对i=123"••”取yXxt）-yXx^）=『y\x）dx=】|）/（“）+y"（G）|Axy,y（兀）一y（£_］）=J：y\x）dx=^［/（%,）+yXx^）］A¥Z,则兀点满足：/（£）=*［y"（£）+y"（兀I）］Ar,+y（£_1）（6）ya）=*［ya）+）‘S-J*+〉"-】）=-{-［/U/）+y"a_］）］Ax,+y（兀_|）+yXx^）}Ax,.+y（xi_］）（7）=-y（兀）心；+-ya_i+y（Vi/+ydi），将⑹、（7）代入（1）,得到ZU）+c心j）+/（%；_!）］Ar,.+y（“）}+〃（兀）1+/（xJAx/++y"（兀T）Axj+y（%/i）Ax.+^^）］=/（x,）从而/“（）一a（兀）丄y（兀T）心,+y（£_］）］—“£•）［：y"（兀T）心；+yXx^）Sxi+y（£_|）］玖戈”--------2--------------------j--------------严一----------------------------------------------⑻1+—a（xi）\xi+—b（£）山「2i=1,2,3,…丿算法：Step1:按（5）式求出/（x0）oStep2:匸0;whilei<n{匸汁1;将y（Vi）、y"（无-）、y"（兀i）代入⑻,输出）产（兀・）;将fg）、y\x.t）代入（6）,输出/（%,.）；将y（s）、y'Cs）、y©）代入歹（兀）=*［：/（兀J+y'Cs）IX+y“（），输出y（兀•）；}3误差分析记/尸心，•二兀-兀i（实际计算中，步长通常是取定的，不随i变化），在矩形算法中，根据（3）式，误差的阶为0（/?）o在梯形算法中，根据（7）式，误差的阶仍为0（h\虽然两种方法误差的阶是相同的，但由于在数值微积分中，没有可以确定的绝对谋差与相对谋差，所以仅凭谋羌的阶还不能完全说明实际误差水平。本文算法对力、y（%i）、#（无t）、丁"...