第八章常微分方程数值解

156第八章常微分方程数值解由常微分方程理论可知，我们只能求一些特殊类型的常微分方程。而实际上许多常微分方程求解非常困难。本章主要讨论一阶常微分方程的初值问题：0,yayyxfdxdybxa（8－1）从理论上讲只要方程中的yxf,连续且关于y满足李普希兹（Lipschitz）条件，即存在常数L，使2121,,yyLyxfyxf则常微分方程存在唯一解)(xyy。微分方程数值解：就是求微分方程的解xy在一系列离散节点bxxxxann110157处的近似值iy（i=1，2，…，n）iiixxh1称为由ix到1ix的步长，通常取为常数h。求数值解，首先将微分方程离散化，常用方法有：（1）用差商代替微商若用向前差商代替微商，即iiiiixyxfxyhxyxy,1（i=1，2，…，n）则得1ixyiiixyxhfxy,即iiiiyxhfyy,1（2）数值积分法利用数值积分法左矩形公式iixyxy1=iixxyxhfdxxyxfii,,1可得同样算法iiiiyxhfyy,1（3）用泰勒（Taylor）公式hxyxyii1iixyhxyiiixyxhfxy,得离散化计算公式iiiiyxhfyy,1158§1欧拉（Euler）方法1.1欧拉方法对一阶微分方程（8—1），等分区间ba,为n份，bxxxxann110，则ihaxinabhni,2,1由以上讨论可知，无论用一阶向前差商，还是用数值积分法左矩形公式，或者用泰勒公式取前两项都可得到同样的离散化计算公式iiiiyxhfyy,1代入初值则得到数值算法：ayyyxhfyyiiii01,（i=1，2，…，n－1）（8－2）称其为欧拉方法。几何上欧拉方法就是用一条折线近似表示曲线xyy。（如图）：159P0PiPi+1O0x1x2xix1ix1.2欧拉方法的误差估计定义1局部截断误差：假设)(iixyy为准确值，用某数值算法计算1iy产生的误差160111iiiyxyR，称为该数值算法的局部截断误差。定义2整体截断误差：准确解)(ixy与数值解iy的误差，iiiyxye。设xy有二阶导数，由泰勒公式有：hxyxyii1=iiiyhxyhxy221=iiiiyhyxhfy221,所以111iiiyxyR=)(212iyh，1iiixx（8－3）当h充分小时，欧拉方法的局部截断误差与h2是同阶无穷小，称其为一阶方法。定义3如果一数值解法的局部截断误差为)(1phO，则称该算法为p阶算法。1.3改进的欧拉方法由微分方程数值解的三种基本构造方法知，若取不同的差商（如向后），不同的161数值积分公式（如梯形公式），以及泰勒公式取前三项、四项等可得不同的算法。如果用梯形公式计算积分：11,,2,1iiiixxxyxfxyxfhdxxyxfii111,,2iiiiiiyxfyxfhyy（8－4）且111iiiyxyR=yh3121（8－5）由于此方程为1iy的隐式方程，不易求解。一般将其与欧拉方法联合使用。可得算法kiiiiikiiiiiyxfyxfhyyyxhfyy111101,,2,（8－6）（k=0，1，2，…；i=1，2，…，n－1）实际计算中，当h比较小时，常取一次迭代后的近似值11iy为1iy，于是有改进的欧拉方162法)]~,(),([2),(~1111iiiiiiiiiiyxfyxfhyyyxhfyy（i=0,1,2,…,n-1）例1用欧拉方法和改进的欧拉方法求微分方程7.0,0,10322xyxydxdy的数值解（取h=0.1）。解：由欧拉方法（8－2），得数值计算公式1iy=iy＋0.1×232iiyx计算结果如表8-1由改进的欧拉方法（8－6），得数值计算公式]~3232[05.0321.0~2112121iiiiiiiiiiyxyxyyyxyy163计算结果如表8-2表8-1xi0.00.10.20.30.40.50.60.7yi1.00001.00691.02081.03911.06281.09231.12691.1643误差0.00000.00370.00770.09930.01200.01510.01890.0222表8-2xi0.00.10.20.30.40.50.60.7yi1.00001.00331.01321.02921.05061.07731.10791.1422误差0.00000.00000.00000.00020.00200.01010.01220.0053例2用欧拉法、改进欧拉法求微分方程数值解（h=0.1）。1)0(yyxyy解由欧拉方法（8－2），得数值计算公式1iy=iy＋0.1iiiyxy由改进的欧拉方法（8－6），得数值计算公式164)]~~()[(05.0)(1.0~11111iiiiiiiiiiiiiyxyyxyyyyxyyy计算结果如下表8-3ixiy（Euler）iyˆ（改进Euler）iiyyˆ00.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.00001.10001.20091.30431.41181.52461.64431.77221.91002.05912.22131.00001.10051.20271.30831.41901.53621.66141.79651.94302.10302.278200.00050.00180.00400.00720.01150.01720.02430.03310.04390.0569§2龙格－库塔（Runge－Kutta）方法2.1泰勒展开法由于欧拉方法为一阶方法，为了提高165算法的阶，有必要讨论更高阶的方法。在泰勒展开式中取更多的项，如取p＋1项可得p阶算法。pipiiiiyphyhyhyy!!221ippiyphR111!1其中kiy）可用复合函数求导法则计算。如p=2时得二阶泰勒方法),(221][2,!2iiyxyxiiiiiiifffhyxhfyyhyhyy2.2龙格－库塔法为了避免计算高阶导数龙格－库塔方法利用yxf,某些点处的值的线性组合构造计算公式，使其按泰勒公式展开后与初值问题解的泰勒展开式比较，有尽可能多的项相同。龙格－库塔法的一般形式为：hyhxfKhyhxfKyxfKKaKaKahyymimimiiimmii,,,222122111（8－7）166下面以二阶龙格－库塔法为例说明龙格－库塔法的构造过程。（8－8）将K2在iiyx,处按泰勒公式展开，则iiiiiyxhfayhKahKayy,122111＋],,,,[2222hOyxfyyxhfyxfxhyxfhaiiiiiiii=322211,,,,hOyxfyyxfyxfxhayxhfaayiiiiiiiii=322221hOxyhaxyhaaxyiii另一方面hxyxyii1=3221hOxyhxyhxyiii所以为使局部截断误差的阶尽可能高，应使1222122111,,hKyhxfKyxfKKaKahyyiiiiii1672112221aaa方程组有无穷多组解，取定参数则得到许多具体的二阶龙格－库塔公式。如取21aa=21，2=1则得龙格－库塔法即改进的欧拉公式。同理可得更高阶的龙格－库塔法。常用的四阶（经典）龙格－库塔法：342312143211,2,22,2,226hKyhxfKkhyhxfKKhyhxfKyxfKKKKKhyyiiiiiiiiii（8－9）例3用二阶龙格－库塔法求例1（5.0,0x，h取0.1）的数值解。解由（8－9）得16823422321221432111.01.03205.005.03205.005.0323222601KyxKKyxKKyxKyxKKKKKyyiiiiiiiiii计算得结果见表8-4表8-4____________________________________________________________________xi0.00.10.20.30.40.5yi1.0000001.0033221.0131591.0291241.0507181.077217误差0.0000000.0000020.0000080.0000110.0000240.000031____________________________________________________________________§3阿达姆斯（Adams）方法将初值问题写成等价形式0011,yxydxxyxfxyxyikixxkii169记xyxfxF,，当k取不同的值以及对xf用不同的插值多项式近似时得到不同的数值算法。3.1阿达姆斯外插公式上式中取k=0，并取321,,,nnnnxxxx为节点，作函数xF的三次插值多项式xp3=30iinxF30ijjjninjnxxxx其插值余项为xR3则xF=xp3＋xR31nxydxxPxynnxxn1积分后得阿达姆斯外插公式nnyy1＋332211,9,37,59,5524nnnnnnnnyxfyxfyxfyxfh（n=3，4，5…）（8－10）其局部截断误差为Rn=nyh557202511,nnnxx（8－11）由于积分区间为1,nnxx，插值区间为nnxx,3，因此称其为阿达姆斯外插公式。1703.2阿达姆斯内插公式同样取k=0，并取211,,,nnnnxxxx为节点，作函数xF的三次插值多项式可得nnyy1＋211519924nnnnffffh（8－12）余项Rnnyh5572019（8－13）由于阿达姆斯内插公式是隐式公式，故用它计算时也需用迭代法。例4用四阶阿达姆斯外插公式求初值问题5.0,0,00xyyxdxdy的数值解，取h=0.1。解由（8－10）得nnyy1＋332211,9,37,59,5524nnnnnnnnyxfyxfyxfyxfh=12.024.09.07.39.55.18241321nyyyynnnn（n=3，4，5）171计算结果见表8-5表8-5______________________________