物理学中常微分方程初值问题的数值解法

1、1 第四章第四章 物理学中常微分方程初值问题的数值解法物理学中常微分方程初值问题的数值解法4.1 物理学中的常微分方程物理学中的常微分方程4.2 常微分方程初值问题的一级、二级欧拉近似法常微分方程初值问题的一级、二级欧拉近似法 4.3 龙格龙格- -库塔法库塔法24.1 物理学中的常微分方程物理学中的常微分方程一、力学中的例子一、力学中的例子l 落体运动落体运动：作用力：作用力重力、阻力重力、阻力牛顿方程：牛顿方程： Kvmgdtdvml阻尼振动阻尼振动：作用力：作用力弹性力、阻尼力弹性力、阻尼力 振动方程：振动方程： vdtdxKvkxdtdvmdtdxKkxdtxdm22一阶常微分方程一阶

2、常微分方程二阶常微分方程二阶常微分方程3)( 020022xdtdxdtxd：阻尼因子：阻尼因子 0：振子固有频率：振子固有频率)( sin200022tFxdtdxdtxd)0( ),0(xv?)( ?)(txtv问题：若已知问题：若已知，求，求l对比阻尼振动的标准形式：对比阻尼振动的标准形式：l受迫振动：受迫振动： 二阶常微分方程二阶常微分方程二阶常微分方程二阶常微分方程4： 弹簧振子弹簧振子：忽略各种阻力和弹簧质量的理想模型：忽略各种阻力和弹簧质量的理想模型 平衡位置：弹簧原长，选为原点平衡位置：弹簧原长，选为原点 ；回复力：；回复力：0,2222xmkdtxdkxdtxdmkmoxf0

3、2022xdtxd单摆：单摆：忽略阻力和摆线质量，摆锤可视为质点，摆角小于忽略阻力和摆线质量，摆锤可视为质点，摆角小于5度度 平衡位置：竖直位置；回复力矩：平衡位置：竖直位置；回复力矩： mglmglsin 由转动定理：由转动定理： 0,22222lgdtdmgldtdml02022dtdmgTloox20mk令kxF令令 lg20二阶齐次线性常微分方程二阶齐次线性常微分方程二阶齐次线性常微分方程二阶齐次线性常微分方程50,2222IcdtdcdtdI由转动定理：由转动定理：Ixy不计阻力和弹簧质量，竖直弹簧振子的运动也是简谐振动不计阻力和弹簧质量，竖直弹簧振子的运动也是简谐振动在平衡位置在平

4、衡位置 , 取为原点取为原点omglk kxxlkmgf )(回复力：回复力：与水平弹簧振子一样也是简谐振动与水平弹簧振子一样也是简谐振动 loxFmg02022dtd20Ic令令扭摆扭摆 ：忽略各种阻力，忽略弹性杆的质量忽略各种阻力，忽略弹性杆的质量 回复力矩回复力矩：c动力学方程：动力学方程：02022xdtxdmk20二阶齐次线性常微分方程二阶齐次线性常微分方程二阶齐次线性常微分方程二阶齐次线性常微分方程6二、电学中的例子二、电学中的例子RC放电电路放电电路上电压：上电压：RRIVR上电压：上电压：CCqVC电路方程：电路方程： 0CRVV（无电源）（无电源） )(dtdqI 0Cqdt

5、dqR（与落体方程相当（与落体方程相当, ,求求 ?)(tq） 0CqRI一阶常微分方程一阶常微分方程RC1、 放电电路放电电路7RLC2、 电磁振荡电磁振荡电路电路电路方程：电路方程： 0CRLVVV（无电源）（无电源）dtdILVL上电压：上电压：LRIVR上电压：上电压：R上电压：上电压：CCqVCIdtdqCqRIdtdIL0022CqdtdqRdtqdL与与阻尼振动阻尼振动方程相当方程相当! ! 问题：问题： ?)( ?)(tItq二阶常微分方程二阶常微分方程8三、三、常微分方程数值解法的原理常微分方程数值解法的原理1 1高阶化为一阶方程高阶化为一阶方程),(),(22ztfdtdz

6、zdtdydtdyytfdtydIvzqxy ( ( ) )两个未知函数两个未知函数)(),(tzty联合求解，联合求解，只需研究一阶方程的解法。只需研究一阶方程的解法。 92. 泰勒级数泰勒级数 设一阶微分方程设一阶微分方程 ),(ytfdtdy的解为的解为 )(ty,则则 ! 3)(! 2)()()(32ttyttytytytty在级数中取若干项，得到近似方法：在级数中取若干项，得到近似方法：一级欧拉法一级欧拉法：取：取2项项 )(2tO （截断误差）（截断误差） 二级欧拉法二级欧拉法：取：取3项项 )(3tO 龙格龙格库塔法库塔法：取：取5项项)(5tO 104.2 常微分方程初值问题的

7、一级、二级欧拉近似法常微分方程初值问题的一级、二级欧拉近似法 一、一、一级欧拉近似法一级欧拉近似法1 1基本公式：基本公式： ),(1tytfytyyyiiiiiititNii ,2 , 1 , 0 求：求： )(ty00)(),(ytyytfdtdy已知已知: :递推公式：递推公式： Nyyyy210112 RC电路电路 0CRVV电路方程：电路方程： 0CqRIdtdqI 求求 ?)(tq0)0( qqtRCqdtdqCqdtdqR无关）与0TRC （ 时间常数）时间常数） (这里 RCqqfytfiiii)(),() )1 (1RCtqtRCqqqiiii12例例1： 300 , 1 ,

8、10 , 1)0(tCRqt取取1， ，求：求： tq 并与精确解并与精确解 RCteqtq0)(比较。比较。 笔算此题：笔算此题： 9 . 010111RCtW9 . 0*1jjjqWqq数值解：数值解： 精确解：精确解： teq*1 . 0)1 (1RCtqtRCqqqiiii133RLC电路电路0022)0( ,)0(0IIqqCqdtdqRdtqdL?)( ?)(tItq求： 数值解数值解： 有关）与 ,( 221IqLCqILRdtqddtdItIqqIdtdqiii)( 1LCqILRftfIIiiiiii而14若若 0)0(I，在，在 CLR/2的条件下，有：的条件下，有： 22

9、000)2/( ,/1 ,cosLRLC其中其中 )cos()(200teqtqtLR解析解解析解：15例例2：已知：已知 1 , 1 , 4 . 0 , 0)0( , 1)0(LCRIq， 200 t, 1 . 0t求：求： )( ),(tItq并与精确解并与精确解 )(tq作比较。作比较。 计算程序：计算程序： real*8 I(0:200),q(0:200),R,C,L,dt write(*,*)input R=?,C=?,L=? read(*,*)R,C,L open(1,file=RCL1.dat) q(0)=1.0 I(0)=0.0 t0=0.0 write(1,*) t0,I(0

10、),q(0) dt=0.1 do 10 K=0,199 f=(1.)*(R*I(K)/L+q(K)/(L*C) I(K+1)=I(K)+f*dt q(K+1)=q(K)+I(K)*dt t=fioat(K+1)*dt10 write(1,*) t, I(K+1), q(K+1) end 221LCqILRdtqddtdItIqqIdtdqiii1, ()iiiiiiqRIIftfILLC 16RLC电路电路一级欧拉近似解一级欧拉近似解)(tq曲线曲线)(tI曲线曲线0246810121416182022-0.8-0.6-0.4-0.20.00.20.40.60.81.0 tq(t) 精确解 一

11、级欧拉近似解17二、二二、二级欧拉近似法级欧拉近似法)(! 2)(321tOttytyyyiiii (1) 求：求： )(ty00)(),(ytyytfdtdy已知已知: :)(!2)(321tOttytyyyiiii (1) )(2121tOtytyyyiiii (2) 将将(2)代入代入(1):18)()(41)(214311tOtOytyytyyyiiiiii (3) 略去略去 )(3tO 项：项： tyytyyyiiiii211)(),(21tOtytfyyiiii(5)注意可用可用一级欧拉法一级欧拉法上式中右边仍有上式中右边仍有 1iy项，但该项有项，但该项有 t相乘，相乘，近似近似

12、,不会降低精度不会降低精度 1iy(4) 2),(),(11tytfytfyiiiii2)(1tyyyiii192),(),(),(1111tytfytfyytytfyyiiiiiiiiii(6) 递推关系递推关系： 11iiiyyy(6)可进一步写为：（由可进一步写为：（由(6)式中的第一式解出式中的第一式解出 ),(iiytf并代入并代入(6)式中的第二式）式中的第二式） 递推公式递推公式：),(21),( 11111tytfyyytytfyyiiiiiiiii20例例3：RC电路：电路：0)0( qqtRCqdtdq无关）与RCqqtfiii),(RCtW1WqtRCqqqiiii*1)

13、(21111tRCqqqqiiii一级欧拉近似结果一级欧拉近似结果二级近似结果：二级近似结果：)*(211Wqqii)1 (212Wqi),(21),(11111tytfyyytytfyyiiiiiiiii21仍以仍以例例1为例，为例， 9 . 010111RCtW，则，则 iiiqWqq*905. 021*2122例例4： RLC电路：电路：0022)0( ,)0(0IIqqCqdtdqRdtqdL求：求： ?)(?)(tItq数值解：数值解： ),( ,( 221ytfdtdyIqLCqILRdtqddtdItIqqIdtdqiii有关）与*),(21)( 11111tItfIIILCqI

14、LRftfIIiiiiiiiiiii而23计算程序计算程序 real*8 I(0:200),q(0:200),R,C,L, dt,I1(0:200) write(*,*)input R=?,C=?,L=? read(*,*)R,C,L open(1,file=RCL2.dat) q(0)=1.0 I(0)=0.0 t0=0.0 write(1,*) t0,I(0),q(0) dt=0.1 do 10 K=0,199 f=-(R*I(K)/L+q(K)/(L*C) I1(K+1)=I(K)+f*dt f1=(1.)*(R*I1(K+1)/L+q(K)/(L*C) I(K+1)=0.5*(I(K)

15、+I1(K+1)+f1*dt) q(K+1)=q(K)+I(K)*dt t=float(K+1)*dt10 write(1,*) t, I(K+1), q(K+1) end一级欧拉近似结果一级欧拉近似结果*),(21)( 11111tItfIIILCqILRftfIIiiiiiiiiiii而二级欧拉近似结果二级欧拉近似结果24U=U0sin t例例5：RL暂态电路暂态电路 一阶一阶RL暂态电路如图所示，其中暂态电路如图所示，其中U0=311V，=314rad，R=10，L=500mH。求开关。求开关K合上后电流合上后电流i(t)在区间在区间0，0.01上的数值解。取上的数值解。取t=0.001

16、。解解：根据电路理论可得初值问题为：根据电路理论可得初值问题为：0)0(sin0iiLRtLUdtdi0)0(20314sin622iitdtdi即：即：25用用一级欧拉近似法一级欧拉近似法公式如下：公式如下：10( , )0.001(622sin31420 )0.980.622sin3140iiiiiiiiiiif t ititiiti用用二级欧拉近似法二级欧拉近似法公式如下：公式如下：11111100.980.622sin314 ()()20.990.311(sin314sin314)0.010iiiiiiiiiiiiiiittiif tif tiittii，26此初值问题的此初值问题的精

17、确解精确解为：为：)cossin()()(220ttLReLRLUtitLRti一级欧拉近似一级欧拉近似二级欧拉近似二级欧拉近似精确解精确解00000.0010.192110.096050.096200.0020.553710.371010.372880.0031.045670.794250.799220.0041.616201.320731.329830.0052.205871.895381.909230.0062.753492.458502.477200.0073.202042.951582.974740.0083.504243.323043.349800.0093.627213.5332

18、33.562350.013.555663.558363.58835计算结果计算结果27则则 x关于时间关于时间 t的微分方程可表示为：的微分方程可表示为： 0222ugmxkdtxd同同例例4类似，将二阶微分方程化成：类似，将二阶微分方程化成： 2mxkugdtdvvdtdx，平方的平方的引力，设引力，设x反比于铁块质心与磁铁中心的距离反比于铁块质心与磁铁中心的距离N块块对地面的压力，开始磁铁与铁块的距离对地面的压力，开始磁铁与铁块的距离为为0.305m，m例例6：如图所示，一块永久磁铁：如图所示，一块永久磁铁A对质量为对质量为的铁块的铁块产生一产生一u，设铁块与地面磨擦系数为，设铁块与地面磨

19、擦系数为k反比系数为反比系数为27200/smk， 1 . 0u2m/s81. 9g， 表示铁表示铁28例例6二级欧拉近似解二级欧拉近似解)(tx曲线曲线)(tv曲线曲线29例例7：导弹追踪问题导弹追踪问题设位于坐标原点的甲舰向位于设位于坐标原点的甲舰向位于x轴上点轴上点A(1, 0)处的乙舰发射处的乙舰发射导弹，导弹头始终对准乙舰导弹，导弹头始终对准乙舰.如果乙舰以最大的速度如果乙舰以最大的速度v0(是常是常数数)沿平行于沿平行于y轴的直线行驶，导弹的速度是轴的直线行驶，导弹的速度是5v0，求导弹运行，求导弹运行的曲线方程，又乙舰行驶多远时，导弹将它击中？的曲线方程，又乙舰行驶多远时，导弹将

20、它击中？解法一（解析法）解法一（解析法）假设导弹在假设导弹在t时刻的位置为时刻的位置为P(x(t), y(t)，乙舰位于，乙舰位于Q(1,v0t)由于导弹头始终对准乙舰，故此时直线由于导弹头始终对准乙舰，故此时直线PQ) 1 ()1 (0yyxtvxytvy10即有：即有：就是导弹的轨迹曲线弧就是导弹的轨迹曲线弧OP在点在点P处的切线处的切线30又根据题意，弧又根据题意，弧OP的长度为的长度为AQ的的5倍倍)2(51002tvdxyx联系联系(1)、(2)式消去式消去t可得：可得：(3) 151)1 (2yyx初值条件为：初值条件为： 0)0(y0)0( y解即为导弹的运行轨迹解即为导弹的运行

21、轨迹 245)1 (125)1 (855654xxy当当 时时 ，即当乙舰航行到点，即当乙舰航行到点 处时被导弹击中。处时被导弹击中。 1x245y)245 , 1 (被击中时间为被击中时间为: 。00245vvyt若若v0=1，则在，则在t=0.21时被击中。时被击中。) 1 ()1 (0yyxtv31解法二（数值方法）解法二（数值方法）令令 ，将方程，将方程(3)化为一阶微分方程。化为一阶微分方程。yz(3) 151)1 (2yyx)1/(1512xzzzy初值条件为：初值条件为： 0)0(y0)0(z由由一阶欧拉近似公式一阶欧拉近似公式可得：可得：0)0()1 (510)0(211zxx

22、zzzyxzyyiiiii32 program main real*8 y(0:200), z(0:200), x0, dx, x open(1,file=DDoula.dat) y(0)=0.0 z(0)=0.0 x0=0.0 write(1,*) x0,y(0),z(0) dx=0.01 do 10 k=0,100 y(k+1)=y(k)+z(k)*dx z(k+1)=z(k)+sqrt(1+z(k)*2)*dx/(5-5*x) x=float(k+1)*dx10 write(1,*) x, y(k+1) end一级欧拉近似一级欧拉近似计算程序如下：计算程序如下：0)0()1 (510)0

23、(211zxxzzzyxzyyiiiii33模拟导弹追踪轨迹：模拟导弹追踪轨迹：试用试用二级欧拉近似二级欧拉近似求解导弹追踪问题？求解导弹追踪问题？0.00.20.40.60.81.00.000.050.100.150.20Y Axis X Axis 34作业作业EX4-1：假设高楼顶上有一物体掉下来，在下落过程中，假设高楼顶上有一物体掉下来，在下落过程中，受到空气的阻力与下落速度的平方成正比，正比系数为受到空气的阻力与下落速度的平方成正比，正比系数为 ，物体质量，物体质量 ，重力加速度为，重力加速度为9.81m/s2，设楼顶处为原点，初速度为设楼顶处为原点，初速度为0，试用，试用一级欧拉近似

24、法、二级一级欧拉近似法、二级欧拉近似法欧拉近似法编程编程求解下落位移、速度随时间的变化关系求解下落位移、速度随时间的变化关系（ ）。）。 mSNK05. 0gm5030 t35tEX4-2：已知弹簧振子已知弹簧振子阻尼振动阻尼振动的方程为：的方程为： )( 020022xdtdxdtxd, 100t时时, 0 . 0 , 0 . 1vx取取0.5秒，试分别取秒，试分别取4 . 0 , 1 . 0用用二级欧拉法二级欧拉法近似画出近似画出 曲线曲线tvtx,)1000(t若是若是受迫振动受迫振动，振动方程为，振动方程为 试取试取 ，求解微分方程并给出，求解微分方程并给出 曲线。曲线。 )( sin

25、20022txdtdxdtxd2 , 1 , 5 . 0tx 36020406080100-6-4-202460=1 tX =0.5 =1 =2受迫振动示意图受迫振动示意图37EX4-3：如右图所示的一单摆，摆长如右图所示的一单摆，摆长 1.016m，重力加速度，重力加速度9.81m/s2，摆，摆角角 所满足的所满足的微分方程可表示为微分方程可表示为求用求用二级欧拉法近似二级欧拉法近似求摆角求摆角 与时间与时间 的曲线关系。的曲线关系。 mll0sin22dtdklmgdtdmlt384.3 龙格龙格- -库塔法库塔法求：求： )(ty00)(),(ytyytfdtdy已知已知: :一、龙格库

26、塔公式一、龙格库塔公式q 类似于欧拉法公式的推导，利用二元函数的类似于欧拉法公式的推导，利用二元函数的Taylor 可以导出一类具有四阶精度的可以导出一类具有四阶精度的龙格龙格-库塔公式库塔公式。 不进行推导不进行推导)(5tO 展开式，使计算公式的局部截断误差取为展开式，使计算公式的局部截断误差取为39常用格式：常用格式：),()2,2()2,2(),(3423121tKyttfKKtyttfKKtyttfKytfKiiiiiiii6)22(43211tKKKKyyii40解：解：精确解：精确解： 3222tteyt1)0() 10( , 12yttyy例例8：取步长：取步长1 . 0h，用

27、，用龙格龙格- -库塔公式库塔公式求解常微分求解常微分方程方程初值问题初值问题)100( ititi输入初值输入初值 10y1 . 0t， , 1),(2tyytf令令 计算得：计算得： 0)1 ,0(),(001fytfK0025. 0) 1 ,05. 0()2,2(1002fKtyttfK002375. 0)2,2(2003KtyttfK0097625. 0),(3004tKyttfK而精确解为：而精确解为： )(1ty000325164. 16/)22() 1 . 0(432101tKKKKyyy000325208. 1计算结果比欧拉方法计算结果比欧拉方法精度要高！精度要高！(见程序）见

28、程序）41 program mainreal*8 y,dydx,xx=0.0d0y=1.0d0dydx=0.0d0h=0.1d0do j=1,10 call rk4(y,dydx,x,h)x=x+hyy=-2*exp(-x)+x*2.-2*x+3write(*,*) real=,x , y,yycall derivs1(x,y,dydx)end do end求下一个点导数求下一个点导数精确解精确解龙格库塔法求解存在龙格库塔法求解存在y中中计算程序计算程序3222tteyt精确解：精确解： 1)0() 10( , 12yttyy42 subroutine rk4(y,dydx,x,h)real*

29、8 y,dydx,yt,dyt,dym,dyk,x,xh,h6,hhhh=h*0.5 ! h6=h/6 !xh=x+hh !2/ t6/ t2tt subroutine derivs1(x,y,dydx) real*8 x,y(1),dydx(1) dydx(1)=-y(1)+x*2.+1 return end 求导数求导数43y=y+h6*(dydx+2*dyt+2*dym+dyk)returnend计算计算y1计算计算K4存在存在dyk计算计算K3存在存在dym计算计算K2存在存在dyt yt=y+hh*dydx call derivs1(xh,yt,dyt) yt=y+hh*dyt ca

30、ll derivs1(xh,yt,dym) yt=y+h*dym call derivs1(x+h,yt,dyk)2001(,)22ttKftyK3002(,)22ttKftyK4003(,)KfttytK 11234(22)6iityyKKKK2001(,)22ttKf tyK3002(,)22ttKf tyK4003(,)Kf tt ytK 44二、常微分方程组的求解二、常微分方程组的求解含有两个未知函数的一阶常微分方程组初值问题含有两个未知函数的一阶常微分方程组初值问题可以写作如下形式：可以写作如下形式： )0(2022122)0(1012111)( ),()( ),(ytyyytfdt

31、dyytyyytfdtdy为了用为了用龙格库塔公式龙格库塔公式求解以上初值问题，求解以上初值问题，方程组写成向量形式，记方程组写成向量形式，记 ， ， 21yyy21yyy),(),(),(212211yytfyytfytf通常将通常将45则则一阶常微分方程组一阶常微分方程组初值问题式可以写为初值问题式可以写为)0(0)(),(ytyytfy)()22()22()(6)22(3)(42)(31)(2)(14321)()1(Ktyt,tfKKty ,ttfKKty ,ttfKy , tfKtKKKKyyiiiiiiiiii龙格库塔公式对龙格库塔公式对一阶微分方程组一阶微分方程组初值问题仍适用初值

32、问题仍适用) 1 (2) 1 (1yy)2(2)2(1yy)(2)(1iiyy) 1(2) 1(1iiyy （向量序列）（向量序列）461)0( , 0)0( , 3 . 0042321212211yytyyyyyy解：此初值问题的解：此初值问题的精确解精确解为：为：)(31)(51tteety)2(31)(52tteety记记 2121123),(yyyytf， 212124),(yyyytf例例9：取步长取步长1 . 0h，用用龙格龙格- -库塔公式库塔公式求解求解常微分方程组初值问题常微分方程组初值问题1) 1 , 0 , 0(),(2) 1 , 0 , 0()0(),0(,( ),(2

33、)0(2)0(1021212101)0(2)0(10111fyytfKfyytfyytfK)()(1iiy , tfK4745. 1)1.05 , 1 . 0 ,05. 0( )2,2,2(4 . 2)1.05 , 1 . 0 ,05. 0( )2,2,2(212)0(211)0(10222112)0(211)0(10121fKhyKhyhtfKfKhyKhyhtfK)22(1)(2Kty ,ttfKii5525. 1)1.0725 ,12. 0 ,05. 0( )2,2,2(505. 2)1.0725 ,0.12 ,05. 0( )2,2,2(222)0(221)0(10232122)0(2

34、21)0(10131fKhyKhyhtfKfKhyKhyhtfK)22(2)(3Kty ,ttfKii48因此有：因此有：247866667. 0)22(6) 1 . 0(41312111)0(11)1(1KKKKhyyy152704167. 1)22(6) 1 . 0(42322212)0(22)1(2KKKKhyyy以上数值解与精确解的误差为以上数值解与精确解的误差为 51046. 9同理可计算出：同理可计算出： 63287176. 0)2 . 0(1)2(1 yy45160267. 1)2 . 0(2)2(2 yy24618565. 1)3 . 0(1)3(1 yy98700407. 1

35、) 3 . 0(2)3(2 yy15725. 2)1.15525 ,2505. 0 , 1 . 0( ),(062. 3)1.15525 ,0.2505 , 1 . 0( ),(232)0(231)0(10242132)0(231)0(10141fhKyhKyhtfKfhKyhKyhtfK)(3)(4Ktyt,tfKii(1)( )1234(22)6iityyKKKK49例例7: 解法三（参数方程）解法三（参数方程） 设时刻设时刻t乙舰的坐标为乙舰的坐标为(X(t)，Y(t)，导弹的坐标为，导弹的坐标为(x(t)，y(t)。3、因乙舰以速度、因乙舰以速度v0沿直线沿直线x=1运动，设运动，设v

36、0=1，则，则w=5，X=1，Y=t1、设导弹速度恒为、设导弹速度恒为w，则，则) 1 ()()(222wdtdydtdx2、由于弹头始终对准乙舰，故导弹的速度平行于、由于弹头始终对准乙舰，故导弹的速度平行于 乙舰与导弹头位置的差向量，乙舰与导弹头位置的差向量， )2(yYxXdtdydtdx即即)3()()()()()()(2222yYyYxXwdtdyxXyYxXwdtdx由由(1),(2)消去消去 可得：可得：500)0(, 0)0()()()1 (5)1 ()()1 (52222yxytytxdtdyxytxdtdx因此导弹运动轨迹的参数方程为：因此导弹运动轨迹的参数方程为：解此参数方

37、程即可得导弹追踪轨迹。解此参数方程即可得导弹追踪轨迹。51三、高阶常微分方程的求解三、高阶常微分方程的求解 m阶微分方程初值问题阶微分方程初值问题 mmmmatyatyatytttyyytfy)( ,.,)( ,)( ),.,(0)1(100010)1()(上式可以转化为一阶微分方程组的初值问题。上式可以转化为一阶微分方程组的初值问题。,1yy,2yy)1(,mmyy 令：令：m阶微分方程初值问题阶微分方程初值问题关于关于 的的 一阶常微分方程组一阶常微分方程组初值问题初值问题 myyy , , ,2152mmmmmmatyatyatyyyytfyyyyyyy)(,)(,)(),(,02021

38、012, 113221例如，考虑二阶微分方程的初值问题：例如，考虑二阶微分方程的初值问题： 100010)( ,)( ),(atyatytttyytfy20210121221)(,)(),(atyatyyytfyyy令令 yy 1yy2一阶常微分方程组一阶常微分方程组初值问题初值问题也可采用降阶方法求解！也可采用降阶方法求解！53例例10：对二：对二阶微分方程阶微分方程初值问题初值问题 6 . 0)0( , 4 . 0)0(10 ,sin222yytteyyyt解：令解：令 ,1yy yy2，化为微分方程组问题：，化为微分方程组问题： 6 . 0)0( , 4 . 0)0(10 22sin21

39、212221yytyyteyyyt, 00t应用应用龙格库塔公式龙格库塔公式，取，取 , 1 . 0h 22sin212yyteft计算可得：计算可得： 6 . 0)0(211 yK4 . 0) 0 (2) 0 (2sin)0 (),0 (,(2102210120yyteyytfKt62. 02)0(12221KhyK3247644757. 0)2)0(,2)0(,2(122111022KhyKhyhtfK)22(1)(2Kty ,ttfKii)()(1iiy , tfK546162382238. 02)0(22231KhyK3152409237. 0)2)0(,2)0(,2(222211032KhyKhyhtfK6315240924. 02)0(32241KhyK2178637298. 0)0(,)0(,(322311042hKyhKyhtfK4617333423. 0)22(6)0() 1 . 0() 1 . 0(4131211111KKKKhyyy6316312421. 0)22(6)0() 1 . 0() 1 . 0(4232221222KKKKhyyy该初值问题的该初值问题的精确解精确解为为)cos2(sin2 . 0)(2t