常微分方程數(shù)值解法縮

1、 第六章常微分方程數(shù)值解法n本章研究常微分方程初值問題 (6.1.1) 的數(shù)值解法. 并且假定f(x, y)滿足解的存在唯一性定理及相當光滑等條件.n初值問題(6.1.1)的精確解記為y(x). d( , ),d( )yf x yaxbxy a建立數(shù)值解法的基本思想n 本章的數(shù)值解法，它不是求(6.1.1)的解y(x)的解析表達式或近似表達式，而是通過某種離散化方法，將連續(xù)變量的初值問題轉(zhuǎn)化為關于離散量的差分方程的初值問題來求一系列離散點上的解值y(xi) 的近似值yi. 利用計算機解微分方程主要使用數(shù)值方法. n取一系列點 x0, x1, , xn , y(x0)=y0 , y(x1)y1

2、, , y(xn)yn , y0 , y1 , , yn , 稱為數(shù)值解用離散化方法建立求y(xn)的近似值yn的遞推格式(差分方程)，求得yn h= xn xn-1稱為步長本章都取定步長.n解初值問題(6.1.1)的數(shù)值解法，其特點都是采取步進式的方法，即求解過程順著節(jié)點排列的次序一步一步向前推進.n這種數(shù)值解法分為兩大類：(1)單步法：若求yn+1 ，只需利用它前一步的 信息yn ，則稱這種方法為單步法。它由 y0出發(fā)，可求得y1 ，y2， y3 (2)多步法：若求yn+1 ，需利用它前面至少兩個點的息，則稱這種方法為多步法.數(shù)值解法研究的主要問題：(1)方法的推導：采用的離散化手段，精度

3、準則.(2)收斂性：差分方程的解是否充分逼近初值問題的解.(3)穩(wěn)定性：初始數(shù)據(jù)、計算過程中每步產(chǎn)生的誤差對以后各步解的影響，這種誤差傳播是否在衰減.n具體的數(shù)值方法還應考慮(1)誤差估計(2)解的起動方法(3)步長如何選取(4)隱式方法的如何計算建立數(shù)值解法的基本途徑n常用的離散化方法 (1)Taylor展開 (2)化導數(shù)為差商 (3)數(shù)值積分6.1 單步法及基本概念Euler折線法n利用Taylor展開法將y(xn+1)在xn處Taylor展開得差分方程 yn+1=yn+hf(xn,yn) （n=0,1,2,）稱此方法為Euler折線法或矩形法1()()(, ()nnnny xy xhf

4、xy x21()()()()()2!nnnnnhy xy xhy xhy xyn 利用化導數(shù)為差商的方法得差分方程 yn+1=yn+hf(xn,yn) （n=0,1,2,） (6.1.2)(,()()()(1nnnnnxyxfxyhxyxy1()()(, ()nnnny xy xhf xy xn利用數(shù)值積分的方法在xn,xn+1上對y(x)=f(x, y (x)積分得用左矩形求積公式計算定積分有 y(xn+1)y(xn)+hf(xn,y(xn)以此得差分方程 yn+1=yn+hf(xn,yn) （n=0,1,2,）1d)(,()()(1nnxxnnxxyxfxyxy梯形法n用不同的近似公式計算

5、定積分的值,就得到解初值問題的不同數(shù)值解法.n用梯形求積公式計算積分得 (6.1.3) 這個方法稱為梯形法. 它是隱格式.n運用它常采用下面的迭代格式111 (,)(,)(0,1,2,)2nnnnnnhyyf xyf xyn(0)1(1)( )111(,) (,)(,)2(0,1,2,;0,1,2,)nnnnkknnnnnnyyhf xyhyyf xyf xykn改進的Euler方法n若梯形法只迭代一次,便得改進的Euler方法 (6.1.4)也可寫成形式 (6.1.4)11(1)(0)11(,)(0,1,2,) (,)(,)2nnnnnnnnnnyyhf xynhyyf xyf xy1121

6、211122(,)(,)nnnnnnyyhkhkkf xykf xh yhk數(shù)值方法精度的衡量準則n定義定義 設yn=y(xn)，則稱Tn=y(xn+1)-yn+1為方法的從xn到xn+1這一步的局部截斷誤差. n定義若差分方程對所有y(x)Mr都精確成立，而對于某個r+1次多項式不能精確成立，則稱這個數(shù)值方法是r階的. 等價定義等價定義 若數(shù)值方法的局部截斷誤差為O(hr+1)，則稱這種方法為r階方法，這里r為非負整數(shù).n方法的階數(shù)越高，逼近效果越好.n由Taylor展開式知 (6.1.5)對于Euler法Tn=y(xn+1)- y(xn)-hf(xn,y(xn)=y(xn+1)-y(xn)

7、-hy(xn) Euler法是1階方法 )(! 2)( )()()(21nnnnnxyhxhyxyhxyxy2()2!nhyn對于梯形法由Taylor展開式知 (6.1.6)(6.1.5)與(6.1.6)相減并利用y(xn+1)=y (xn)+O(h)梯形法法是2階方法.11()()()()(0,1,2,)2nnnnnhTy xy xy xy xn2111()()()()()2!nnnnnhy xy xhy xhy xy x311()()()()()2nnnnnhTy xy xy xy xO hn對于改進Euler法，設yn=y(xn)，利用Taylor展開式nK1=f(xn,y(xn)=y(

8、xn) 211(, ()(, ()(, ()(, ()(, ()(, ()()(, ()()()nnnnnnnnnnnnnnnnnKf xh y xhkf xy xhf xy xhkf xy xxyf xy xhf xy xy xf xy xxyy xhy x將其代入中Tn有改進Euler法是2階方法.31121()()()()2nnnTy xy xkkO hRunge-Kutta方法nTaylor展開式n理論上講，只要解y(x)充分光滑，通過保留Taylor展開式的若干項就可得到任意階的近似公式，但計算y(x)的各階導數(shù)很麻煩?？砷g接利用這種思想 )(,()(,()(,(! 2)(,()()

9、(! 2)( )()()(221nnnnnnnnnnnnnnxyxfyxyxfxyxfxhxyxhfxyxyhxhyxyhxyxynEuler法也可寫成形式n其局部截斷誤差為O(h2) ，是一階方法.每步計算f的值一次111(,)nnnnyyk hkf xynEuler預報-校正公式也可寫成形式n局部截斷誤差為O(h3)，是二階方法每步計算f的值二次n可以通過增加計算f的值的次數(shù)，提高公式的階數(shù)（精度） 1121211122(,)(,)nnnnnnyyhkhkkf xykf xh yk hn以f在不同點上的函數(shù)值的線性組合來代替yn+1 yn，其中有一些可待定選取的待定參數(shù)，通過Taylor展

10、開確定這些待定參數(shù)使建立的數(shù)值方法按要求達到一定的階數(shù)，這種思想就是Runge-Kutta方法的思想nRunge-Kutta法的一般形式其中Ri,ai,bij都是常數(shù)，a0=0，b1j=0， ( j=1,2,i-1)1111(,),1,2, .snniiiiininijjjyyhR Kkf xa h yhb Kisn二階Runge-Kutta法的一般形式 其中R1, R2, a, b為待定常數(shù) 其局部截斷誤差為O(h3)，是二階方法每步計算f的值二次11122121(,)(,)nnnnnnyyRhkR hkKf xyKf xah ybhkn設yn=y(xn) K1=f(xn,yn)=f(xn,

11、y(xn)把K2中f在(xn,y(xn)處Taylor展開21212(, ()(, ()(, ()(, ()()(, ()(, ()(, ()(, ()()nnnnnnnnnnnnnnnnKf xah y xbhKf xy xahf xy xbhKf xy xO hxyf xy xh af xy xbf xy xf xy xO hxyn再將K1,K2代入yn+1中，n將其與y(xn+1)泰勒展開式比較，要使y(xn+1)-yn+1=O(h3)，含h0, h1, h2的項相同即有11122122322()() (, ()(, ()(, ()(, ()()nnnnnnnnnnnyyRhKR hKy

12、 xh RRf xy xh aRf xy xbR f xy xf xy xO hxy個未知數(shù)，個方程滿足條件的解不止一組取R1=R2=1/2,a=b=1,就是改進Euler法；取R1=0,R2=1,a=1/2,b=1/2的中點方法.122211212RRaRbRn四階Runge-Kutta法的一般形式其中有13個待定常數(shù)局部截斷誤差為O(h)，是四階方法每步計算f的值四次1112233441222113331132244411422433(,)(,)(,)(,)nnnnnnnnnnyyRhKR hKR hKR hKKf xyKf xa h yb hkKf xa h yb hkb hkKf xa

13、 h yb hkb hkb hkn設yn=y(xn)n把K2,K3,K4中f在(xn,y(xn)處泰勒展開后，將K1,K2,K3,K4代入yn+1中，再將yn+1按h的冪重新整理后與y(xn+1)泰勒展開式比較，要使y(xn+1)-yn+1=O(h5)，含h0, h1, h2 , h3, h4的項相同從而確定各參數(shù)n 24181)(121)(61)(4131211443322433422443323243234222433222433422433224343332324243232224433224342414323132124321RbbababaRaRbaababaRRbababaRRba

14、RaRaRaRaRaRaRaRaRabbbabbabaRRRR13個未知數(shù)，11個方程滿足條件的解不止一組最常用的是n標準四階Runge-Kutta公式 (6.1.7)112341122343(22)6(,)(,)22(,)22(,)nnnnnnnnnnhyyKKKKKf xyhKhKf xyhKhKf xyKf xh yhKn單步法的一般形式為 yn+1=yn+h(xn,yn,h) (6.1.8)其中多元函數(shù)依賴于f.單步法的收斂性nen= y(xn) yn稱為整體截斷誤差. 收斂性就是討論當x=xn固定且h=(xn-x0)/n0時 en 0的問題.n定義定義 如果對f(x, y)滿足解的存

15、在唯一性條件的初值問題(6.1.1)，差分方程 (6.1.8)的解對每個確定xa,b滿足則稱單步法(6.1.8)是收斂的.0lim( ),nhnyy xnhxan定義定義 若增量函數(shù)(x,y,h)使 (x,y,0)=f(x,y)成立，則稱單步法(6.1.8)與(6.1.1) 相容.n滿足相容條件的方法至少是1階的.n定理定理 設增量函數(shù)(x,y,h)在區(qū)域axb, -y, 0hh0中連續(xù)，并且對變量y滿足Lipschitz條件, 即(x,y1,h)-(x,y2,h)Ly1-y2在這個前提下單步法(6.1.8)收斂的充分必要條件是相容性條件成立.單步法的絕對穩(wěn)定性n當步長取定后，計算中的誤差隨著

16、步數(shù)的增加會不會積累到超出我們許可的范圍，這就是穩(wěn)定性問題.n單步法(6.1.8)應用于模型方程 y=y (0) 設得到的解為yn+1=R(h)yn當0時,實驗方程的精確解y(x)=e(x-a)按模遞減的,這要求 (6.1.8)的解yn是遞減的，誤差也是遞減的，即要求滿足R(h) 1n定義定義 單步法(6.1.8)應用于模型方程 y=y(0) ，若得到的解yn+1=R(h)yn滿 足R(h) 1, 則稱單步法(6.1.8)是絕對穩(wěn)定 的.使R(h) 1成立的 =h所在區(qū)間稱為 絕對穩(wěn)定區(qū)間.n例 標準四階Runge-Kutta公式 由R(h) 0時,實驗方程的精確解y(x)=e(x-a)按模遞

17、增的; (6.4.1)的解當n時無界；誤差也是無界的.這種情況下，如果誤差相對于真解是小的，就說方法是相對穩(wěn)定的. 這時要求 ri(h)r0(h) ，i=1,p.n當0時,實驗方程的精確解y(x)=e(x-a)按模遞減的; 要求(6.4.1)的解當n時遞減到0，誤差也是遞減的.這種情況下，就需要ri(h)1, i=0,1,p來保證.這時就稱方法是絕對穩(wěn)定的. 方法(6.2.1)的穩(wěn)定性取決于特征多項式(r; h)零點 的性質(zhì)，故又稱(r; h)為(6.2.1)的穩(wěn)定多項式. 穩(wěn)定性的定義 記=hn定義 設方法(6.2.1)是收斂的. ri() (i=0,1,p)是穩(wěn)定多項式(r; )=0的根,

18、 r0()是形如r0(h)=1+h+O(h2)的根.(1)若對任意的,R有 ri()r0() ，i=1,p.且當ri()=r0() 時, ri()是單根,則稱方法(6.2.1)在,上是相對穩(wěn)定的， ,稱為此方法的相對穩(wěn)定區(qū)間. (2)若對任意的(,)R有 ri()0時,相對穩(wěn)定性重要； 當0時,絕對穩(wěn)定重要，才有意義.n定義 若一個方法的絕對穩(wěn)定區(qū)間是 (-,0), 則稱該方法是A穩(wěn)定的.n(1)收斂的方法的相對穩(wěn)定區(qū)間不空。因為方法穩(wěn)定的前提是收斂，故滿足根條件和相容條件,0,.n(2)只有0的情況討論絕對收斂性才有意義.n(3)從誤差分析的角度看，絕對穩(wěn)定的方法是理想的.n(3)相對穩(wěn)定區(qū)

19、間和絕對穩(wěn)定區(qū)間越大越好，可使盡可能大的一類微分方程是穩(wěn)定的.例 討論梯形方法的相對和絕對穩(wěn)定區(qū)間.n解 梯形方法為對于試驗方程y=y ，梯形方法成為即其穩(wěn)定多形式為它只有一個根，記為且r0(0)=r0=1( ; )(1)(1)22hhr hr1(1)(1)22nnhhyy11()2nnnnhyyyy11()2nnnnhyyyy012( )12hr hh 由于由于其只有一個根，所以梯形方法的相對穩(wěn)定區(qū)間是(-，+) 當 0時, r0() 1 當 0時, r0()1 這樣，梯形方法的絕對穩(wěn)定區(qū)間是(-，0).即梯形方法是A穩(wěn)定的.例 討論Simpson方法的相對和絕對穩(wěn)定區(qū)間nSimpson方法

20、為利用該格式求解試驗方程y=y 得其穩(wěn)定多形式為它有兩個根24( ; )(1)(1)333hhr hrhr1111(4)3nnnnnhyyyyy20021133( )(0,( )1),113hhr hhr hh2121133( )113hhr hh114(1)(1)333nnnhhyhyy考察當 0 時，上式小于等于1；當 0時,當 0時,n因此，Simpson方法不存在絕對穩(wěn)定區(qū)間.2021211333|( )|1111133hhhr hhh2112211333| ( )|1111133hhhr hhhn步長的選取:穩(wěn)定性與步長有關，步長的選取一定要保證方法是穩(wěn)定的，即h屬于穩(wěn)定區(qū)間.n收斂

21、性是反映遞推公式本身的整體截斷誤差對計算結(jié)果的影響；穩(wěn)定性反映某一計算步驟中出現(xiàn)的誤差對計算結(jié)果的影。6.5 預測-校正方法線性多步法(6.2.1)當b-10時，是隱格式估計出yn+1的一個初值y(0)n+1，利用迭代y(j)n+1是yn+1的第j次近似值.當hb-11時，此迭代收斂. 1110()pnnin iin iiyb hfa yhb f(1)( )11110(,)()pjjnnnin iin iiyb hf xya yhb fn用顯式方法來作預測值y(0)n+1 ，用隱式方法迭代校正一次得y(1)n+1值，這種顯式與隱式聯(lián)合使用構成的方法稱為預報-校正法.作為預報的顯式公式稱為預測式

22、，用于校正的隱式公式稱為校正式.n例如，改進的Euler方法.11(1)(0)11(,)(0,1,2,) (,)(,)2nnnnnnnnnnyyhf xynhyyf xyf xyn預測公式與校正公式選取同階方法，可使截斷誤差用預測值和校正值表示，再用截斷誤差來修正，得到提高精度的方法.n4階顯式Adams方法 局部截斷誤差為n4階隱式Adams方法 局部截斷誤差為 1112919524nnnnnnhyyffff1123555937924nnnnnnhyyffff5(5)11251()()720nnnnTy xyh yx5(5)1119()()720nnnnTy xyh yx n4階Adams預報-校正方法稱為4階Adams預報-校正格式(PECE模式).1(0)123(0)(0)111(1)(0)1112(1)(1)111:555937924:(,):919524:(,)nnnnnnnnnnnnnnnnnnhP yyffffEff xyhC yyffffEff xy(0)5(5)11(1)5(5)11251()()72019()()720nnnnnny xyh yxy xyh yx (0)(0