物理學(xué)中常微分方程初值問題的數(shù)值解法

1、1 第四章第四章 物理學(xué)中常微分方程初值問題的數(shù)值解法物理學(xué)中常微分方程初值問題的數(shù)值解法4.1 物理學(xué)中的常微分方程物理學(xué)中的常微分方程4.2 常微分方程初值問題的一級、二級歐拉近似法常微分方程初值問題的一級、二級歐拉近似法 4.3 龍格龍格- -庫塔法庫塔法24.1 物理學(xué)中的常微分方程物理學(xué)中的常微分方程一、力學(xué)中的例子一、力學(xué)中的例子l 落體運動落體運動：作用力：作用力重力、阻力重力、阻力牛頓方程：牛頓方程： Kvmgdtdvml阻尼振動阻尼振動：作用力：作用力彈性力、阻尼力彈性力、阻尼力 振動方程：振動方程： vdtdxKvkxdtdvmdtdxKkxdtxdm22一階常微分方程一階

2、常微分方程二階常微分方程二階常微分方程3)( 020022xdtdxdtxd：阻尼因子：阻尼因子 0：振子固有頻率：振子固有頻率)( sin200022tFxdtdxdtxd)0( ),0(xv?)( ?)(txtv問題：若已知問題：若已知，求，求l對比阻尼振動的標(biāo)準(zhǔn)形式：對比阻尼振動的標(biāo)準(zhǔn)形式：l受迫振動：受迫振動： 二階常微分方程二階常微分方程二階常微分方程二階常微分方程4： 彈簧振子彈簧振子：忽略各種阻力和彈簧質(zhì)量的理想模型：忽略各種阻力和彈簧質(zhì)量的理想模型 平衡位置：彈簧原長，選為原點平衡位置：彈簧原長，選為原點 ；回復(fù)力：；回復(fù)力：0,2222xmkdtxdkxdtxdmkmoxf0

3、2022xdtxd單擺：單擺：忽略阻力和擺線質(zhì)量，擺錘可視為質(zhì)點，擺角小于忽略阻力和擺線質(zhì)量，擺錘可視為質(zhì)點，擺角小于5度度 平衡位置：豎直位置；回復(fù)力矩：平衡位置：豎直位置；回復(fù)力矩： mglmglsin 由轉(zhuǎn)動定理：由轉(zhuǎn)動定理： 0,22222lgdtdmgldtdml02022dtdmgTloox20mk令kxF令令 lg20二階齊次線性常微分方程二階齊次線性常微分方程二階齊次線性常微分方程二階齊次線性常微分方程50,2222IcdtdcdtdI由轉(zhuǎn)動定理：由轉(zhuǎn)動定理：Ixy不計阻力和彈簧質(zhì)量，豎直彈簧振子的運動也是簡諧振動不計阻力和彈簧質(zhì)量，豎直彈簧振子的運動也是簡諧振動在平衡位置在平

4、衡位置 , 取為原點取為原點omglk kxxlkmgf )(回復(fù)力：回復(fù)力：與水平彈簧振子一樣也是簡諧振動與水平彈簧振子一樣也是簡諧振動 loxFmg02022dtd20Ic令令扭擺扭擺 ：忽略各種阻力，忽略彈性桿的質(zhì)量忽略各種阻力，忽略彈性桿的質(zhì)量 回復(fù)力矩回復(fù)力矩：c動力學(xué)方程：動力學(xué)方程：02022xdtxdmk20二階齊次線性常微分方程二階齊次線性常微分方程二階齊次線性常微分方程二階齊次線性常微分方程6二、電學(xué)中的例子二、電學(xué)中的例子RC放電電路放電電路上電壓：上電壓：RRIVR上電壓：上電壓：CCqVC電路方程：電路方程： 0CRVV（無電源）（無電源） )(dtdqI 0Cqdt

5、dqR（與落體方程相當(dāng)（與落體方程相當(dāng), ,求求 ?)(tq） 0CqRI一階常微分方程一階常微分方程RC1、 放電電路放電電路7RLC2、 電磁振蕩電磁振蕩電路電路電路方程：電路方程： 0CRLVVV（無電源）（無電源）dtdILVL上電壓：上電壓：LRIVR上電壓：上電壓：R上電壓：上電壓：CCqVCIdtdqCqRIdtdIL0022CqdtdqRdtqdL與與阻尼振動阻尼振動方程相當(dāng)方程相當(dāng)! ! 問題：問題： ?)( ?)(tItq二階常微分方程二階常微分方程8三、三、常微分方程數(shù)值解法的原理常微分方程數(shù)值解法的原理1 1高階化為一階方程高階化為一階方程),(),(22ztfdtdz

6、zdtdydtdyytfdtydIvzqxy ( ( ) )兩個未知函數(shù)兩個未知函數(shù))(),(tzty聯(lián)合求解，聯(lián)合求解，只需研究一階方程的解法。只需研究一階方程的解法。 92. 泰勒級數(shù)泰勒級數(shù) 設(shè)一階微分方程設(shè)一階微分方程 ),(ytfdtdy的解為的解為 )(ty,則則 ! 3)(! 2)()()(32ttyttytytytty在級數(shù)中取若干項，得到近似方法：在級數(shù)中取若干項，得到近似方法：一級歐拉法一級歐拉法：?。喝?項項 )(2tO （截斷誤差）（截斷誤差） 二級歐拉法二級歐拉法：?。喝?項項 )(3tO 龍格龍格庫塔法庫塔法：?。喝?項項)(5tO 104.2 常微分方程初值問題的

7、一級、二級歐拉近似法常微分方程初值問題的一級、二級歐拉近似法 一、一、一級歐拉近似法一級歐拉近似法1 1基本公式：基本公式： ),(1tytfytyyyiiiiiititNii ,2 , 1 , 0 求：求： )(ty00)(),(ytyytfdtdy已知已知: :遞推公式：遞推公式： Nyyyy210112 RC電路電路 0CRVV電路方程：電路方程： 0CqRIdtdqI 求求 ?)(tq0)0( qqtRCqdtdqCqdtdqR無關(guān)）與0TRC （ 時間常數(shù)）時間常數(shù)） (這里 RCqqfytfiiii)(),() )1 (1RCtqtRCqqqiiii12例例1： 300 , 1 ,

8、10 , 1)0(tCRqt取取1， ，求：求： tq 并與精確解并與精確解 RCteqtq0)(比較。比較。 筆算此題：筆算此題： 9 . 010111RCtW9 . 0*1jjjqWqq數(shù)值解：數(shù)值解： 精確解：精確解： teq*1 . 0)1 (1RCtqtRCqqqiiii133RLC電路電路0022)0( ,)0(0IIqqCqdtdqRdtqdL?)( ?)(tItq求： 數(shù)值解數(shù)值解： 有關(guān)）與 ,( 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) 遞推關(guān)系遞推關(guān)系： 11iiiyyy(6)可進(jìn)一步寫為：（由可進(jìn)一步寫為：（由(6)式中的第一式解出式中的第一式解出 ),(iiytf并代入并代入(6)式中的第二式）式中的第二式） 遞推公式遞推公式：),(21),( 11111tytfyyytytfyyiiiiiiiii20例例3：RC電路：電路：0)0( qqtRCqdtdq無關(guān)）與RCqqtfiii),(RCtW1WqtRCqqqiiii*1)

13、(21111tRCqqqqiiii一級歐拉近似結(jié)果一級歐拉近似結(jié)果二級近似結(jié)果：二級近似結(jié)果：)*(211Wqqii)1 (212Wqi),(21),(11111tytfyyytytfyyiiiiiiiii21仍以仍以例例1為例，為例， 9 . 010111RCtW，則，則 iiiqWqq*905. 021*2122例例4： RLC電路：電路：0022)0( ,)0(0IIqqCqdtdqRdtqdL求：求： ?)(?)(tItq數(shù)值解：數(shù)值解： ),( ,( 221ytfdtdyIqLCqILRdtqddtdItIqqIdtdqiii有關(guān)）與*),(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一級歐拉近似結(jié)果一級歐拉近似結(jié)果*),(21)( 11111tItfIIILCqILRftfIIiiiiiiiiiii而二級歐拉近似結(jié)果二級歐拉近似結(jié)果24U=U0sin t例例5：RL暫態(tài)電路暫態(tài)電路 一階一階RL暫態(tài)電路如圖所示，其中暫態(tài)電路如圖所示，其中U0=311V，=314rad，R=10，L=500mH。求開關(guān)。求開關(guān)K合上后電流合上后電流i(t)在區(qū)間在區(qū)間0，0.01上的數(shù)值解。取上的數(shù)值解。取t=0.001

16、。解解：根據(jù)電路理論可得初值問題為：根據(jù)電路理論可得初值問題為：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計算結(jié)果計算結(jié)果27則則 x關(guān)于時間關(guān)于時間 t的微分方程可表示為：的微分方程可表示為： 0222ugmxkdtxd同同例例4類似，將二階微分方程化成：類似，將二階微分方程化成： 2mxkugdtdvvdtdx，平方的平方的引力，設(shè)引力，設(shè)x反比于鐵塊質(zhì)心與磁鐵中心的距離反比于鐵塊質(zhì)心與磁鐵中心的距離N塊塊對地面的壓力，開始磁鐵與鐵塊的距離對地面的壓力，開始磁鐵與鐵塊的距離為為0.305m，m例例6：如圖所示，一塊永久磁鐵：如圖所示，一塊永久磁鐵A對質(zhì)量為對質(zhì)量為的鐵塊的鐵塊產(chǎn)生一產(chǎn)生一u，設(shè)鐵塊與地面磨擦系數(shù)為，設(shè)鐵塊與地面磨

19、擦系數(shù)為k反比系數(shù)為反比系數(shù)為27200/smk， 1 . 0u2m/s81. 9g， 表示鐵表示鐵28例例6二級歐拉近似解二級歐拉近似解)(tx曲線曲線)(tv曲線曲線29例例7：導(dǎo)彈追蹤問題導(dǎo)彈追蹤問題設(shè)位于坐標(biāo)原點的甲艦向位于設(shè)位于坐標(biāo)原點的甲艦向位于x軸上點軸上點A(1, 0)處的乙艦發(fā)射處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對準(zhǔn)乙艦導(dǎo)彈，導(dǎo)彈頭始終對準(zhǔn)乙艦.如果乙艦以最大的速度如果乙艦以最大的速度v0(是常是常數(shù)數(shù))沿平行于沿平行于y軸的直線行駛，導(dǎo)彈的速度是軸的直線行駛，導(dǎo)彈的速度是5v0，求導(dǎo)彈運行，求導(dǎo)彈運行的曲線方程，又乙艦行駛多遠(yuǎn)時，導(dǎo)彈將它擊中？的曲線方程，又乙艦行駛多遠(yuǎn)時，導(dǎo)彈將

20、它擊中？解法一（解析法）解法一（解析法）假設(shè)導(dǎo)彈在假設(shè)導(dǎo)彈在t時刻的位置為時刻的位置為P(x(t), y(t)，乙艦位于，乙艦位于Q(1,v0t)由于導(dǎo)彈頭始終對準(zhǔn)乙艦，故此時直線由于導(dǎo)彈頭始終對準(zhǔn)乙艦，故此時直線PQ) 1 ()1 (0yyxtvxytvy10即有：即有：就是導(dǎo)彈的軌跡曲線弧就是導(dǎo)彈的軌跡曲線弧OP在點在點P處的切線處的切線30又根據(jù)題意，弧又根據(jù)題意，弧OP的長度為的長度為AQ的的5倍倍)2(51002tvdxyx聯(lián)系聯(lián)系(1)、(2)式消去式消去t可得：可得：(3) 151)1 (2yyx初值條件為：初值條件為： 0)0(y0)0( y解即為導(dǎo)彈的運行軌跡解即為導(dǎo)彈的運行

21、軌跡 245)1 (125)1 (855654xxy當(dāng)當(dāng) 時時 ，即當(dāng)乙艦航行到點，即當(dāng)乙艦航行到點 處時被導(dǎo)彈擊中。處時被導(dǎo)彈擊中。 1x245y)245 , 1 (被擊中時間為被擊中時間為: 。00245vvyt若若v0=1，則在，則在t=0.21時被擊中。時被擊中。) 1 ()1 (0yyxtv31解法二（數(shù)值方法）解法二（數(shù)值方法）令令 ，將方程，將方程(3)化為一階微分方程?；癁橐浑A微分方程。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模擬導(dǎo)彈追蹤軌跡：模擬導(dǎo)彈追蹤軌跡：試用試用二級歐拉近似二級歐拉近似求解導(dǎo)彈追蹤問題？求解導(dǎo)彈追蹤問題？0.00.20.40.60.81.00.000.050.100.150.20Y Axis X Axis 34作業(yè)作業(yè)EX4-1：假設(shè)高樓頂上有一物體掉下來，在下落過程中，假設(shè)高樓頂上有一物體掉下來，在下落過程中，受到空氣的阻力與下落速度的平方成正比，正比系數(shù)為受到空氣的阻力與下落速度的平方成正比，正比系數(shù)為 ，物體質(zhì)量，物體質(zhì)量 ，重力加速度為，重力加速度為9.81m/s2，設(shè)樓頂處為原點，初速度為設(shè)樓頂處為原點，初速度為0，試用，試用一級歐拉近似

24、法、二級一級歐拉近似法、二級歐拉近似法歐拉近似法編程編程求解下落位移、速度隨時間的變化關(guān)系求解下落位移、速度隨時間的變化關(guān)系（ ）。）。 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，擺，擺角角 所滿足的所滿足的微分方程可表示為微分方程可表示為求用求用二級歐拉法近似二級歐拉法近似求擺角求擺角 與時間與時間 的曲線關(guān)系。的曲線關(guān)系。 mll0sin22dtdklmgdtdmlt384.3 龍格龍格- -庫塔法庫塔法求：求： )(ty00)(),(ytyytfdtdy已知已知: :一、龍格庫

26、塔公式一、龍格庫塔公式q 類似于歐拉法公式的推導(dǎo)，利用二元函數(shù)的類似于歐拉法公式的推導(dǎo)，利用二元函數(shù)的Taylor 可以導(dǎo)出一類具有四階精度的可以導(dǎo)出一類具有四階精度的龍格龍格-庫塔公式庫塔公式。 不進(jìn)行推導(dǎo)不進(jìn)行推導(dǎo))(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計算結(jié)果比歐拉方法計算結(jié)果比歐拉方法精度要高！精度要高！(見程序）見

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求下一個點導(dǎo)數(shù)求下一個點導(dǎo)數(shù)精確解精確解龍格庫塔法求解存在龍格庫塔法求解存在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 求導(dǎo)數(shù)求導(dǎo)數(shù)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二、常微分方程組的求解二、常微分方程組的求解含有兩個未知函數(shù)的一階常微分方程組初值問題含有兩個未知函數(shù)的一階常微分方程組初值問題可以寫作如下形式：可以寫作如下形式： )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以上數(shù)值解與精確解的誤差為以上數(shù)值解與精確解的誤差為 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: 解法三（參數(shù)方程）解法三（參數(shù)方程） 設(shè)時刻設(shè)時刻t乙艦的坐標(biāo)為乙艦的坐標(biāo)為(X(t)，Y(t)，導(dǎo)彈的坐標(biāo)為，導(dǎo)彈的坐標(biāo)為(x(t)，y(t)。3、因乙艦以速度、因乙艦以速度v0沿直線沿直線x=1運動，設(shè)運動，設(shè)v

36、0=1，則，則w=5，X=1，Y=t1、設(shè)導(dǎo)彈速度恒為、設(shè)導(dǎo)彈速度恒為w，則，則) 1 ()()(222wdtdydtdx2、由于彈頭始終對準(zhǔn)乙艦，故導(dǎo)彈的速度平行于、由于彈頭始終對準(zhǔn)乙艦，故導(dǎo)彈的速度平行于 乙艦與導(dǎo)彈頭位置的差向量，乙艦與導(dǎo)彈頭位置的差向量， )2(yYxXdtdydtdx即即)3()()()()()()(2222yYyYxXwdtdyxXyYxXwdtdx由由(1),(2)消去消去 可得：可得：500)0(, 0)0()()()1 (5)1 ()()1 (52222yxytytxdtdyxytxdtdx因此導(dǎo)彈運動軌跡的參數(shù)方程為：因此導(dǎo)彈運動軌跡的參數(shù)方程為：解此參數(shù)方

37、程即可得導(dǎo)彈追蹤軌跡。解此參數(shù)方程即可得導(dǎo)彈追蹤軌跡。51三、高階常微分方程的求解三、高階常微分方程的求解 m階微分方程初值問題階微分方程初值問題 mmmmatyatyatytttyyytfy)( ,.,)( ,)( ),.,(0)1(100010)1()(上式可以轉(zhuǎn)化為一階微分方程組的初值問題。上式可以轉(zhuǎn)化為一階微分方程組的初值問題。,1yy,2yy)1(,mmyy 令：令：m階微分方程初值問題階微分方程初值問題關(guān)于關(guān)于 的的 一階常微分方程組一階常微分方程組初值問題初值問題 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應(yīng)用應(yīng)用龍格庫塔公式龍格庫塔公式，取，取 , 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