1、一階微分方程求解講義3.3 一階微分方程的求解 000)(),()( ytytttyfty一階微分方程的求解可歸結(jié)為在給定初始條件下，一階微分方程的求解可歸結(jié)為在給定初始條件下，求微分方程的初值問題求微分方程的初值問題 基本思想：基本思想： 在初值問題存在唯一解的時間區(qū)間內(nèi)，在若干個時間在初值問題存在唯一解的時間區(qū)間內(nèi)，在若干個時間離散點上，用差分方程代替微分方程，然后逐點求解差分離散點上，用差分方程代替微分方程，然后逐點求解差分方程，得到各時間離散點方程，得到各時間離散點 、 處的函數(shù)處的函數(shù) 近似近似值值 、 1t2tnt1y2yny)(ty 當兩相鄰離散點之間的間隔較小時，用一階差商當兩

2、相鄰離散點之間的間隔較小時，用一階差商取代一階導數(shù)取代一階導數(shù) 一一. .前向歐拉法前向歐拉法)( )()(11kkkkktytttyty htytytykkk)( )()(1 hyyykkk1 kktth 1令步長令步長 ，則，則其近似值為：其近似值為： 近似解的誤差首先是由差商差商近似代替微商微商引起的，這種近似代替所產(chǎn)生的誤差稱為截斷誤差。還有一種誤差稱為舍入誤差，這種誤差是由于計算時數(shù)值舍入引起的。前向歐拉法的幾何意義： y(t) y1 y2 y3 y(t3) y(t2) y(t1) y0 y(t0) h h h t0 0 t1 t2 t3 在任一步長內(nèi)，用一段直在任一步長內(nèi)，用一段直

3、線代替函數(shù)線代替函數(shù) 的曲線，的曲線，此直線段的斜率等于該函此直線段的斜率等于該函數(shù)在該步長起點的斜率。數(shù)在該步長起點的斜率。 )(ty),(1kkkkythfyy 歐拉法的幾何意義：過點A0（t0，y0），A1（t1，y1），A n（t n，y n ），斜率分別為f（t0，y0），f（t1，y1），f（tn，y n）所連接的一條折線，所以歐拉法亦稱為歐拉折線法歐拉折線法。例例1. 1. 應用前向歐拉法解初值問題應用前向歐拉法解初值問題0)1(, 21 ,22 ytetytyt)2(21ntnnnnnetythyy nnhttyyn1 . 010)1(00 取步長取步長h=0.1，并把計算結(jié)果

4、與精確解比較，并把計算結(jié)果與精確解比較解：據(jù)前向歐拉法解：據(jù)前向歐拉法又又有：有：271828183. 0)2(0200001 tetythyy684755578. 0)2(1211112 tetythyy【思路】【思路】 用歐拉法求解常微分方程的初值問題時，首先熟練掌握歐拉公式的用歐拉法求解常微分方程的初值問題時，首先熟練掌握歐拉公式的一般形式，根據(jù)具體題目寫出找出歐拉公式的迭代式，并根據(jù)初始條件和所一般形式，根據(jù)具體題目寫出找出歐拉公式的迭代式，并根據(jù)初始條件和所給步長進行迭代求解。給步長進行迭代求解。微分方程微分方程 是一階線性微分方程，是一階線性微分方程，可求出其通解：可求出其通解：

5、tetyty22 C)(2tety )-(2eetyt )-()(2eettyntnn 則方程的解為：則方程的解為： 從而有：從而有： 345919876. 0)(1211 eetyt86642536. 0)(2222 eetyt帶入初值帶入初值 可得可得 0) 1 ( ye C一階非齊次線性微分方程一階非齊次線性微分方程計算結(jié)果列表（計算結(jié)果列表（ 為前向歐拉法計算近似值，為前向歐拉法計算近似值， 為精確值）為精確值） n 01.000011.10.2718281830.3459198760.07401969321.20.6847555780.86664253603

6、1.2769783441.6072150790.33023673541.42.0935476882.6203595520.52681186451.53.1874451223.9676662950.78022117361.64.6208178465.7209615271.10014368171.76.4663963787.9638734791.497477101ntny)(ntynnyty )(ny)(nty正正分析：分析：當步長不是很小時，前向歐拉法的精度不當步長不是很小時，前向歐拉法的精度不是很高。步長取定后，步數(shù)越多，誤差越是很高。步長取定后，步數(shù)越多，誤差越大。大。二、后向歐拉法)( )

7、()(11 kkktyhtyty歐歐拉拉隱隱式式公公式式 hyyykkk11用一階差商近似代替用一階差商近似代替 在一個步長終點的一階導在一個步長終點的一階導數(shù)，則原微分方程化為：數(shù)，則原微分方程化為：),()( ttyfty )(ty00)(yty 對于給定初始條件對于給定初始條件的微分方程的微分方程其近似值：其近似值： 在任一步長內(nèi)，用一段直線在任一步長內(nèi)，用一段直線代替函數(shù)代替函數(shù) 的曲線，此直的曲線，此直線段的斜率等于該函數(shù)在該線段的斜率等于該函數(shù)在該步長終點的斜率。步長終點的斜率。 后向歐拉法的幾何意義： y(t) y1 y2 y3 y(t3) y(t2) y(t1) y0 y(t0

8、) h h h t0 0 t1 t2 t3 t y(t) ),(111 kkkkythfyy精確值精確值近似值近似值注：后向歐拉法的兩種處理方式注：后向歐拉法的兩種處理方式 前向前向Euler法為顯式，后向法為顯式，后向Euler法法為隱式為隱式須解出須解出yk+1. 可用迭代法可用迭代法yk+1 (n+1) = yk + hf (tk+1，yk+1(n)n = 0，1，2， 解得解得yk+1 ，其中其中yk+1(0) = yk + hf (tk，yk).（結(jié)合前（結(jié)合前向歐拉法，預報）向歐拉法，預報）例例2. 2. 應用后向歐拉法解初值問題應用后向歐拉法解初值問題0)1(,21 ,22 yt

9、etytyt12112111121:)2(11即即ntnnntnnnnnthehtyyetythyynn nnhttyyn1 . 010)1(00 取步長取步長h=0.1，并把計算結(jié)果與精確解比較，并把計算結(jié)果與精確解比較解：據(jù)后向歐拉法解：據(jù)后向歐拉法又計算結(jié)果列表（計算結(jié)果列表（ 為后向歐拉法計算近似值，為后向歐拉法計算近似值， 為精確值）為精確值） n 01.000011.10.4442827750.345919876-0.09836289921.21.1068555350.866642536-0.24021299931.32.0409606121.607215079-0.4337455

10、3341.43.3084097732.620359552-0.68805022151.54.9809113233.967666295-1.01324502861.67.1415858565.720961527-1.42062432971.79.8866975397.963873479-1.922824060ntny)(ntynnyty )(ny)(nty負負三. 梯形法及其預估-矯正法 000)(),()( ytyttytfty)( )( 21)()(11 kkkktytyhtyty)( )( 2)()(11 kkkktytyhtyty用一階差商近似地代替函數(shù)在一個步長起點和終點的用一階差商近

11、似地代替函數(shù)在一個步長起點和終點的一階導數(shù)的平均值一階導數(shù)的平均值 梯形公式梯形公式（歐拉中點公式）（歐拉中點公式）近似值：近似值：)(211 kkkkyyhyy改進歐拉法改進歐拉法顯然，梯形公式是隱式法，一般求顯然，梯形公式是隱式法，一般求 需要解方程，需要解方程，常采用迭代法，初值由顯式的歐拉公式給出：常采用迭代法，初值由顯式的歐拉公式給出：1 ky),()0(1kkkkythfyy ),(),(2)0(11)1(1 kkkkkkytfytfhyy然后將然后將 替代梯形公式等式右邊出現(xiàn)的替代梯形公式等式右邊出現(xiàn)的)0(1 ky1 ky),(),(2)n(11)1(1kkkkknkytfyt

12、fhyy當步長當步長h足夠小，且由前向歐拉法計算的已是較好足夠小，且由前向歐拉法計算的已是較好的近似，則迭代一、二次即可的近似，則迭代一、二次即可, 2 , 1 , 0n預報預報校正校正迭代次數(shù)迭代次數(shù)l幾何意義幾何意義lEuler法法折線法折線法l改進改進Euler法法平均斜率折線法平均斜率折線法例例3. 應用梯形預估應用梯形預估-矯正法解初值問題矯正法解初值問題0)1(,21 ,22 ytetytyt取步長取步長h=0.1，并把計算結(jié)果與精確解比較，并把計算結(jié)果與精確解比較kktkkkkketythyyy201 . 0)2 . 01(1 ）（解：據(jù)前向歐拉法解：據(jù)前向歐拉法梯形預估梯形預估

13、- -矯正矯正2205. 0),(),(21121)0(12)0(111）（）（ kkktkktkkkkkkkkkketytetytyytfytfhyy計算結(jié)果列表（計算結(jié)果列表（ 為梯形預估為梯形預估- -矯正法計算矯正法計算近似值，近似值， 為精確值）為精確值）ky)(nty k 01.000011.10.3423777890.3459198760.00354208721.20.8583145370.8666425360.00832799931.31.5927496431.6072150790.01446543641.42.5982982392.6203595520.02206131351

14、.53.9364441143.9676662950.03122218161.65.6789071035.7209615270.04205442471.77.9092092167.9638734790.054664263ktky)(ktykkyty)(function T Y=Trapezia_reckon(odefun,ab,ya,M)% odefun: 微分方程微分方程 a、b：計算區(qū)間：計算區(qū)間% ya：初值：初值 y(a) M：等分數(shù)目：等分數(shù)目% T: 離散的時間變量離散的時間變量 Y梯形公式的預估校正法解梯形公式的預估校正法解h=(ab(2)-ab(1)/M; 步長步長T=zeros

15、(1,M+1);Y=zeros(1,M+1);T=ab(1):h:ab(2); Y(1)=ya;for j=1:M k1=feval(odefun,T(j),Y(j); k2=feval(odefun,T(j+1),Y(j)+h*k1); Y(j+1)=Y(j)+(h/2)*(k1+k2);end )(2)0(11 jjjjyyhyyFunction y=euler_3_3_2(t,x)y=2/t*x+t2*exp(t ）T Y=Trapezia_reckon ( euler_3_3_2,1 2,0,10) 求解器求解器求解問題求解問題特點特點說明說明ode45非剛性非剛性一步算法；一步算法；

16、4，5階階Runge-Kutta算法算法大部分場合的首選算法大部分場合的首選算法ode23非剛性非剛性一步算法；一步算法；2，3階階Runge-Kutta算法算法使用于精度較低的情形使用于精度較低的情形ode113非剛性非剛性多步法；變階次的多步法；變階次的Adams -Bashforth- Moulton 算法算法計算時間比計算時間比ode45短短ode23t剛性剛性采用梯形算法采用梯形算法適合中度剛性問題的求解適合中度剛性問題的求解ode15s剛性剛性多步法；采用了數(shù)值差多步法；采用了數(shù)值差分算法分算法若若ode45失效時，可嘗試失效時，可嘗試使用使用;ode23s剛性剛性一步法；一步法；2階階Rosebrock算法算法當精度較低時，計算時間當精度較低時，計算時間比比ode15s短短ode23tb剛性剛性隱式隱式Runge-Kutta算法算法當精度較低時，計算時間當精度較低時，計算時間比比ode15s短短不同求解器的特點不同求解器的特點 在用常微分方程描述一個電路的暫態(tài)過程時，往往又包含著多個變化速度相差十分懸殊的子過程，這樣一類過程就認為具有“剛性(stiff)”，描述這類過程的微分方程稱為“剛性問