數(shù)值分析71常微分方程初值問題數(shù)值解引言

第7章常微分方程初值問題數(shù)值解法7.1引言7.2歐拉方法與改進(jìn)的歐拉方法7.3龍格-庫塔方法7.1引言本章將考察一階方程的初值問題

我們知道，只要f(x,y)適當(dāng)光滑—譬如關(guān)于y滿足利普希茨(Lipschitz)條件理論上就可以保證初值問題的解y＝f(x)存在并且唯一.

雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來求解一些特殊類型的方程，實際問題中歸結(jié)出來的微分方程主要靠數(shù)值解法.

所謂數(shù)值解法,就是尋求解y(x)在一系列離散節(jié)點上的近似值y1,y2,,yn,yn+1,.相鄰兩個節(jié)點的間距hn=xn+1-xn稱為步長.假定hi=h(i=1,2,)為定數(shù),這時節(jié)點為xi=x0+ih(i=0,1,2,)(等距節(jié)點).

初值問題的數(shù)值解法有個基本特點:都采取“步進(jìn)式”，即求解過程順著節(jié)點排列的次序一步一步地向前推進(jìn).

首先，要對微分方程離散化，建立求解數(shù)值解的遞推公式.一類是計算yn+1時只用到前一點的值yn，稱為單步法.另一類是用到y(tǒng)n+1前面k點的值yn,yn-1,,yn-k+1，稱為k步法.其次，要研究公式的局部截斷誤差和階，數(shù)值解yn與精確解y(xn)的誤差估計及收斂性，還有遞推公式的計算穩(wěn)定性等問題.7.2歐拉方法與改進(jìn)的歐拉方法7.2.1歐拉法與后退歐拉法在xoy平面上，微分方程(1.1)式的解y=f(x)稱作它的積分曲線，積分曲線上一點(x,y)的切線斜率等于函數(shù)f(x,y)的值.如果按f(x,y)在xoy平面上建立一個方向場，那么，積分曲線上每一點的切線方向均與方向場在該點的方向相一致.基于上述幾何解釋，我們從初始點P0(x0,y0)出發(fā),先依方向場在該點的方向推進(jìn)到x=x1上一點P1，然后再從P1點依方向場在該點的方向推進(jìn)到x=x2

上一點P2,循環(huán)前進(jìn)做出一條折線P0P1P2.

一般地，設(shè)已做出該折線的頂點Pn,過Pn(xn,yn)依方向場的方向再推進(jìn)到Pn+1(xn+1,yn+1)，顯然兩個頂點Pn,Pn+1的坐標(biāo)有關(guān)系這就是著名的(顯式)歐拉(Euler)公式.若初值y0已知，則依公式(2.1)可逐次逐步算出各點數(shù)值解.即

例1

用歐拉公式求解初值問題

解取步長h=0.1，歐拉公式的具體形式為其中xn=nh=0.1n(n=0,1,,10),已知y0=1,由此式可得依次計算下去，部分計算結(jié)果見下表.與準(zhǔn)確解相比，可看出歐拉公式的計算結(jié)果精度很差.

xn

歐拉公式數(shù)值解yn準(zhǔn)確解y(xn)

誤差0.20.40.60.81.01.1918181.3582131.5089661.6497831.7847701.1832161.3416411.4832401.6124521.7320510.0086020.0165720.0257260.0373310.052719

歐拉公式具有明顯的幾何意義,就是用折線近似代替方程的解曲線，因而常稱公式(2.1)為歐拉折線法.

還可以通過幾何直觀來考察歐拉方法的精度.假設(shè)yn=y(xn),即頂點Pn落在積分曲線y=y(x)上，那么，按歐拉方法做出的折線PnPn+1便是y=y(x)過點Pn的切線.從圖形上看,這樣定出的頂點Pn+1偏離了原來的積分曲線，可見歐拉方法是相當(dāng)粗糙的.

為了分析計算公式的精度，通?？捎锰├照归_將y(xn+1)在xn處展開，則有在yn=y(xn)的前提下，f(xn,yn)=f(xn,y(xn))=y(xn).于是可得歐拉法(2.1)的公式誤差為稱為此方法的局部截斷誤差.

如果對方程(1.1)從xn到xn+1積分，得右端積分用左矩形公式hf(xn,y(xn))近似，再以yn代替y(xn)，yn+1代替y(xn+1)也得到歐拉公式(2.1)，局部截斷誤差也是(2.3).稱為(隱式)后退的歐拉公式.

如果右端積分用右矩形公式hf(xn+1,y(xn+1))近似，則得到另一個公式

后退的歐拉公式與歐拉公式有著本質(zhì)的區(qū)別,后者是關(guān)于yn+1的一個直接計算公式，這類公式稱作是顯式的；前者公式的右端含有未知的yn+1，它實際上是關(guān)于yn+1的一個函數(shù)方程,這類方程稱作是隱式的.

顯式與隱式兩類方法各有特點，考慮到數(shù)值穩(wěn)定性等其他因素，人們有時需要選用隱式方法，但使用顯式算法遠(yuǎn)比隱式方便.注：隱式方程通常用迭代法求解，而迭代過程的實質(zhì)是逐步顯式化.

設(shè)用歐拉公式給出迭代初值

，用它代入(2.5)式的右端，使之轉(zhuǎn)化為顯式，直接計算得然后再用代入(2.5)式，又有如此反復(fù)進(jìn)行，得如果f(x,y)對y滿足Lipschitz條件(1.3).由(2.6)減(2.5)得由此可知，只要hL<1，迭代法(2.6)就收斂到解.7.2.2梯形方法

為得到比歐拉法精度高的計算公式，在等式(2.4)

右端積分用梯形求積公式近似,并用yn代替y(xn),yn+1代替y(xn+1)，則得稱為矩形方法.

矩形方法是隱式單步法，用迭代法求解，同后退的歐拉方法一樣，仍用歐拉法提供迭代初值，則矩形迭代公式為為了分析迭代過程的收斂性,將(2.7)與(2.8)相減,得于是有使得則當(dāng)k→∞時有,這說明迭代過程(2.8)是收斂的.7.2.3改進(jìn)的歐拉公式我們看到，梯形方法雖然提高了精度，但其算法復(fù)雜，在應(yīng)用迭代公式(2.8)進(jìn)行實際計算時，每迭代一次，都要重新計算函數(shù)f(x,y

)的值，而迭代又要反復(fù)進(jìn)行若干次，計算量很大，而且往往難以預(yù)測.為了控制計算量，通常只迭代一兩次就轉(zhuǎn)入下一步的計算，這就簡化了算法.具體地說，我們先用歐拉公式求得一個初步的近似值，稱之為預(yù)測值，此預(yù)測值的精度可能很差，再用梯形公式(2.7)將它校正一次，即按(2.8)式迭代一次，這個結(jié)果稱之為校正值.這樣建立的預(yù)測—校正系統(tǒng)通常稱為改進(jìn)的歐拉公式:或表為下列平均化形式(2.9)預(yù)測校正例2

用改進(jìn)的歐拉法解例1中的初值問題(2.2).

解仍取步長h=0.1,改進(jìn)的歐拉公式為部分計算結(jié)果見下表

xnyn

誤差

xnyn

誤差00.20.41.1.1840961.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.05818同例1中的歐拉法的計算結(jié)果比較，明顯改善了精度.

例3

(兩種方法的精度比較)用歐拉公式和改進(jìn)的歐拉公式解下述初值問題，并與其準(zhǔn)確解y=e-x+x進(jìn)行比較.解取步長h=0.1，xk=kh(k=0,1,,6).用兩種方法進(jìn)行計算對應(yīng)結(jié)果及絕對誤差見下表

xn

歐拉公式

改進(jìn)歐拉公式

yn

誤差

yn

誤差0.00.10.20.30.40.50.61.0000001.0000001.0100001.0290001.0561001.0904901.13144104.8×10-38.7×10-31.2×10-21.4×10-21.6×10-21.8×10-211.0050001.0190251.0412181.0708021.1070761.14940401.6×10-42.9×10-44.0×10-44.8×10-45.5×10-45.9×10-47.2.3單步法的局部截斷誤差與階

初值問題(1.1),(1.2)的單步法可用一般形式表示為其中多元函數(shù)與f(x,y

)有關(guān)，當(dāng)含有yn+1時，方法是隱式的，若不含yn+1則為顯式方法，所以顯式單步法可表示為(x,y,h)稱為增量函數(shù)，例如對歐拉法(2.1)有它的局部截斷誤差已由(2.3)給出,對一般顯式單步法則可如下定義.

定義1

設(shè)y(x)是初值問題(1.1),(1.2)的準(zhǔn)確解,稱為顯式單步法(2.11)的局部截斷誤差.

Tn+1之所以稱為局部的，是假設(shè)在xn前各步?jīng)]有誤差.當(dāng)yn=y(xn)時，計算一步，則有

所以，局部截斷誤差可理解為用方法(2.11)計算一步的誤差，也即公式(2.11)中用準(zhǔn)確解y(x)代替數(shù)值解產(chǎn)生的公式誤差.根據(jù)定義,顯然歐拉法的局部截斷誤差為即為(2.3)的結(jié)果.這里稱為局部截斷誤差主項.顯然Tn+1=O(h2).一般情形的定義如下：定義2

設(shè)y(x)是初值問題的準(zhǔn)確解，若存