《計算方法》課件1第8章

第八章常微分方程初值問題的數(shù)值解法8.1引言8.2歐拉方法及其改進8.3龍格—庫塔方法8.4線性多步法8.5算法的穩(wěn)定性及收斂性8.6一階常微分方程組與高階方程8.7解微分方程的波形松弛方法8.8微分方程邊值問題的數(shù)值方法習題88.1引言

在高等數(shù)學中，對于常微分方程的求解，給出了一些典型方程求解析解的基本方法，如可分離變量法、常系數(shù)齊次線性方程的解法、常系數(shù)非齊次線性方程的解法等。但能求

解的常微分方程仍然是有限的，大多數(shù)的常微分方程是不可能給出解析解的。另外，在許多實際問題中，并不需要方程解的表達式，而僅僅需要獲得解在若干點上的近似值即可。因此，研究常微分方程的數(shù)值解法很有必要。在本章里，將向大家著重介紹一階常微分方程初值問題：

(8.1)

在區(qū)間［a,b］上的數(shù)值解法。

下面介紹一些與數(shù)值解法有關的內(nèi)容。

定理8.1如果函數(shù)f(x，y)在帶形區(qū)域D={(x,y)|a≤x≤b，

y∈R}上有定義且連續(xù)，并且滿足李普希茲條件，即當

(x,

y1)，(x，y2)∈D時，有

|f(x，y1)－f(x，y2)|≤L|y1－y2|

其中，L(0<L<+∞)稱為李普希茲常數(shù)(它與x,y無關)，則對任何(x0,y0)∈D，初值問題((8.1)式)在［a,b］上存在唯一連續(xù)可微解y=y(x)。

定理8.2如果函數(shù)f(x，y)在帶形區(qū)域D={(x，y)|a

≤x≤b，y∈R}上滿足李普希茲條件，則初值問題((8.1)式)的數(shù)值解是穩(wěn)定的。

(證略。)

所謂常微分方程初值問題的數(shù)值解法，就是要算出精確解y=y(x)在區(qū)間［a,b］上的一系列離散節(jié)點a≤x0<x1<x2<…

<xn<…≤b處的函數(shù)值y(x0),y(x1),…，y(xn)的近似值y0，y1，…，yn，…。相鄰兩個節(jié)點的間距Δxi=xi－xi－1稱為步長，步長可以相等，也可以不等。在本章中，一般取等步長并用字母h表示，并且記f(xn，y(xn))=y′(xn)。求初值問題的數(shù)值解法采用的是“步進式”，即求解過程順著節(jié)點排列的次序一步一步地向前推進，算出yn后再算yn+1。這種數(shù)值方法有單步法和多步法之分。單步法是在計算yn+1

時，僅利用yn；多步法是在計算yn+1時不僅要用到y(tǒng)n，還要用到y(tǒng)n－1，yn－2,…，一般k步法要用到y(tǒng)n,yn－1,…，

yn－k+1

。單步法的代表方法是龍格—庫塔方法；多步法的代表是亞當斯方法。

單步法和多步法又有顯式和隱式之分。顯式單步法形如yn+1=yn+hφ(xn,yn,h)；隱式單步法形如yn+1=yn+hφ(xn,yn,

yn+1,h)。對多步法來說，顯式與隱式方法與單步法相似。

8.2歐拉方法及其改進

8.2.1歐拉公式

歐拉(Euler)方法是解初值問題((8.1)式)最簡單的數(shù)值方法。推導歐拉公式有多種方法，譬如利用差商、數(shù)值積分、數(shù)值微分、泰勒展開法等。下面我們以數(shù)值積分為例進行推導。對(8.1)式的第一式y(tǒng)′=f(x,y)兩端在區(qū)間［xn,xn+1］上進行積分，得

即

(8.2)選擇不同的計算方法計算積分項,

就會得到一系列不同形式的歐拉公式。

1.向前歐拉法

用左矩形公式計算積分項，得

代入(8.2)式中，并用yn近似代替式中y(xn),即可得到向前歐拉公式:

(8.3)

這是一種顯式形式的方程,因此也稱為顯式歐拉公式。

2.向后歐拉法

用右矩形公式計算積分項，得

代入(8.2)式中，并用yn+1近似代替式中y(xn+1)，即可得向后歐拉公式:

yn+1=yn+hf(xn+1,yn+1)

(8.4)

這是一種隱式形式的方程,因此也稱為隱式歐拉公式。由于數(shù)值積分的矩形方法精度很低，因此歐拉公式比較粗糙。

歐拉方法的幾何意義十分清楚，如圖8-1所示。(8.1)式的解曲線y=y(x)過P0(x0,y0)點，從P0出發(fā)以f(x0,y0)為斜率做直線段，與x=x1相交于P1(x1，y1)，顯然有y1=y0+hf(x0,y0)。同理，再從P1出發(fā)，以f(x1,y1)為斜率做直線段，與x=x2相交于P2(x2,y2)。依次類推，可得一條折線P0P1…Pn…，作為解曲線y=y(x)的近似曲線，故歐拉方法又稱為歐拉折線法。圖8-1

3.梯形方法

為了提高精度，用梯形公式計算積分項,

得

代入(8.2)式中，并用yn近似代替式中y(xn)，yn+1近似代替式中y(xn+1)，即可得到梯形公式:

(8.5)由于數(shù)值積分的梯形公式要比矩形公式精度高，因此梯形公式((8.5)式)要比歐拉公式((8.3)式)精度高。梯形法是一種隱式單步法，從n=0開始，每步都要解關于yn+1的一個方程。一般來說，這是一個非線性方程，因此要用迭代法求解。

4.改進歐拉方法

顯式歐拉公式計算工作量小，但精度低。梯形公式雖提高了精度，但為隱式公式，需用迭代法求解，計算工作量大。綜合歐拉公式和梯形公式便可得到改進的歐拉公式。

具體地說，就是先用歐拉公式求出一個初步的近似值

y［0］n+1，稱之為預估值。預估值的精度可能很差，再用梯形公式將它校正一次，即迭代一次，求得的yn+1稱為校正值，由這種預測—校正方法得到的公式稱為改進的歐拉公式：

(8.6)其中，第一式稱為預估算式；第二式稱為校正算式。有時為了計算方便，將(8.6)式改寫為按照上述思想，如果預估時用其他公式(譬如中矩形公式、辛普森公式等)，則可以得到其他形式的預估—校正式。例如，若采用中矩形公式，則有下述公式：8.2.2單步法的局部截斷誤差和階

定義8.1稱

en=y(xn)－yn

為某一數(shù)值方法在點xn處的整體截斷誤差；稱

Rn,h=y(xn+1)－y(xn)－h(huán)φ(xn,y(xn),h)

為顯式單步法式y(tǒng)n+1=yn+hφ(xn,yn,h)在xn+1處的局部截斷誤差。其中，y(x)為初值問題的精確解。

Rn，h之所以稱為局部的，是因為如果假設y(xn)=yn，即第n步以及以前各步都沒有誤差，則由式y(tǒng)n+1=yn+hφ(xn,yn,h)所得的yn+1與y(xn+1)之差為

即在假定的y(xn)=yn條件下，Rn，h=y(xn+1)－yn+1，這就是

Rn,h稱為局部的含義。

定義8.2若數(shù)值方法的局部截斷誤差為O(hp+1)，則稱這種方法為p階的。

通常p越大，h越小，則截斷誤差越小，數(shù)值方法越精確。

可以證明，歐拉方法是一階方法，改進的歐拉方法以及梯形方法是二階方法。下面僅給出歐拉方法的證明。

設y(xn)=yn，把y(xn+1)在xn處展開成泰勒公式，即則由歐拉公式可得

兩式相減得歐拉公式的局部截斷誤差為若y(x)在［a,b］上充分光滑，且令,

則

故歐拉方法是一階方法。

例8.1應用向前歐拉公式求初值問題:

取步長h=0.1，將計算結(jié)果與精確解y=x+e-x對照。

解將區(qū)間［0,1］進行10等分，h=0.1，xn=nh(n=0.1,

…，10)。

向前歐拉公式為

數(shù)值解yn與精確解y(xn)及誤差列于表8.1。

例8.2用改進歐拉法求初值問題：

取步長h=0.1，將計算結(jié)果與精確解y=e-x比較。(0≤x≤1)

解將區(qū)間［0,1］進行10等分，步長h=0.1，節(jié)點xn=nh(n=0.1,…,10)。

應用改進的歐拉公式，得

由此，y1=0.905,y2=0.9052,…，y10=0.90510。

數(shù)值解yn與精確解y(xn)及誤差列于表8.2。

例8.3有初值問題：

試證明：

(1)用梯形公式求得的數(shù)值解為

(2)當步長h→0時，yn收斂于精確解y=e-x。(0≤x≤1)

證明

(1)應用梯形公式:

整理上式，得由此公式遞推可得

于是

(2)設區(qū)間是等距劃分的，對于任意給定的節(jié)點x=xn=nh，步長。顯然當h→0的同時，n→∞，

由此

證畢。8.3龍格—庫塔方法

本節(jié)中將向大家介紹求解初值問題的一類高精度的單步法——龍格—庫塔(Runge-Kutta)方法。8.3.1龍格—庫塔方法的基本思想

我們首先從歐拉公式及改進的歐拉公式著手進行分析。

歐拉公式可改寫為

用它計算yn+1需要計算一次f(x，y)的值。若設yn=y(xn)，則yn+1的表達式與y(xn+1)在xn處的泰勒展開式的前兩項完全相同，即局部截斷誤差為O(h2)。改進的歐拉公式又可改寫為

用它計算yn+1需要計算兩次f(x，y)的值。若設yn=y(xn)，則yn+1的表達式與y(xn+1)在xn處的泰勒展開式的前三項完全相同，即局部截斷誤差為O(h3)。這兩組公式在形式上有一個共同點：都是用f(x，y)在某些點上值的線性組合得出y(xn+1)的近似值yn+1，而且增加計算f(x，y)的次數(shù)，可提高截斷誤差的階。

因此可考慮用函數(shù)f(x，y)在若干點上的函數(shù)值的線性組合來構造近似公式，構造時要求近似公式在(xn，yn)處的泰勒展開式與解y(x)在xn處的泰勒展開式的前面幾項重合，從而獲得達到一定精度的數(shù)值計算公式。這就是龍格—庫塔的基本思想。8.3.2龍格—庫塔方法的推導

按照上述思想，龍格—庫塔方法的一般形式設定為

(8.7)

其中，ωi，αi，βij為待定常數(shù)。從上式可以看到用它計算yn+1需要計算r次f(x，y)的值，因此(8.7)式稱為r階R-K方法。

下面我們來了解幾種常用的R-K方法。

1.二階龍格—庫塔方法

當r=2時，龍格—庫塔格式為

適當選取參數(shù)ω1，ω2，α2，β21的值，使得在yn=y(xn)的假設下，局部截斷誤差為y(xn+1)－yn+1=O(h3)。為此把K1，K2代入第一式y(tǒng)n+1=yn+ω1K1+ω2K2中，得然后在(xn，yn)處作泰勒展開，可得

(8.8)

這里記fn=f(xn,yn)。而微分方程y′=f(x,y)的精確解y=y(x)在點xn處的泰勒展開式為

(8.9)把(8.8)式與(8.9)式進行比較，為使y(xn+1)-yn+1=O(h3),令h,h2項系數(shù)相等，則有下列方程組:

此方程有無窮多個解，從而有無窮多個二階龍格—庫塔方法。按此方程組解出的每一組解，對應的二階龍格—庫塔方法都是二階的。常用的二階龍格—庫塔方法有

(1)當ω1=ω2=

,α2=β21=1時，有

這正好是改進的歐拉公式。

(2)當ω1=0,ω2=1,α2=β21=時，有

這種方法稱為中間點法。

2.三階龍格—庫塔方法

此時r=3，一般格式為把K1，K2，K3代入yn+1的表達式中，再在(xn，yn)處作泰勒展開，然后與y(xn+1)在xn點處的泰勒展開式作比較，并且要使y(xn+1)－yn+1=O(h4)?？傻孟铝蟹匠探M：

此方程有無窮多個解，因此有無窮多個三階龍格—庫塔方法。常用的三階龍格—庫塔方法有

此方法稱為三階庫塔方法。

3.四階龍格—庫塔方法

在實際應用中，最常用的龍格—庫塔方法是四階龍格—庫塔方法。類似于二、三階龍格—庫塔方法的推導，可得一個含有13個未知量、11個方程的方程組。由于推導復雜，這里從略。此處只介紹最常用的一種四階龍格—庫塔方法：

此方法稱為經(jīng)典四階龍格—庫塔方法。

例8.4取步長h=0.2，用經(jīng)典四階R-K公式求解初值問題:

解分析：按照龍格—庫塔法的解題步驟應先算出各Ki(i=1，2，3，4)，再代入

中計算。已知f(x，y)=2xy，x0=0，y0=1，h=0.2，由四階龍格—庫塔公式可得

代入yn+1=yn+

(K1+2K2+2K3+K4)(其中n=0,1,2,…，5)求解。本例方程的解為

，數(shù)值解yn與精確解y(xn)的對照表如表8.3所示。龍格—庫塔方法的推導基于泰勒展開方法，因而它要求所求的解具有較好的光滑性。如果解的光滑性差，那么，使用四階龍格—庫塔方法求得的數(shù)值解，其精度可能反而不如改進的歐拉方法。在實際計算時，應當針對問題的具體特點選擇合適的算法。

8.4線性多步法

8.4.1線性多步法的基本思想

單步法在計算yn+1時，只用到前一步的信息yn，而沒有用到前面幾步計算得到的信息。此時要提高精度，需重新計算多個點處的函數(shù)值，如R-K方法、計算量較大。若能充分

利用第n+1步前面的多步信息來計算yn+1，便可獲得更高精度的算法，這就是多步法的基本思想。多步法中最常用的是線性多步法。如果記y(xn)的近似值為yn，xn=x0+nh，并記fn=f(xn，yn)，則k步線性多步法的一般形式為

(8.10)

其中，αj，βj(j=0，1，…，k)都是實常數(shù)，且αk≠0，|α0|+|β0|≠0。一般可設αk=1，此時(8.10)式可寫為

(8.11)如果βk≠0，(8.10)式就是隱式方程；如果βk=0，它就是顯式方程。如果k=1，(8.10)式就是一種單步法。比如：

k=1,

α0=－1,β0=1,β1=0時，就是歐拉方法；k=1，

α0=－1，β0=，β1=時，就是梯形方法。

下面介紹線性多步法的階和局部截斷誤差的定義。

對于一般的線性多步法((8.11)式)，定義算子R為

(8.12)其中,y(x)為區(qū)間［a,b］上任意一個連續(xù)可微函數(shù)。若y(x)充分可微，將y(x+jh)及y′(x+jh)按泰勒式展開有

把它們代入(8.12)式可得

(8.13)其中,(8.14)給出如下定義：

定義8.3若在(8.13)式中

c0=c1=…=cp=0，cp+1≠0

則稱算子R［y(x)；h］對應的線性多步法((8.10)式)是p階

方法。

定義8.4設y(x)是初值問題的解，(8.10)式是一種線性多步法，則

稱為是(8.10)式在xn+k處的局部截斷誤差。Rn+k按h展開的首項稱為主局部截斷誤差。

若(8.10)式是一種p階線性多步法，則由(8.13)式可知

Rn+k=cp+1hp+1y(p+1)(xn)+O(hp+2)

因此，主局部截斷誤差就是cp+1hp+1y(p+1)(xn)，而cp+1稱為主局部截斷誤差系數(shù)。8.4.2線性多步法的構造

構造線性多步法有多種方法，常用的是數(shù)值積分法和泰勒展開法(待定系數(shù)法)。

1.數(shù)值積分法

初值問題可以寫成與之等價的積分方程形式設{xi}為等距節(jié)點，xi+1=xi+h。拉格朗日插值基函數(shù)li(t)為

(i=0，1，2，…，p)

用拉格朗日插值多項式來近似代替(8.14)式中的被積函數(shù)，得到近似計算公式：

若用yn-i代替y(xn－i)，用fn－i表示f(xn－i,yn－i)，則可得線性多步法顯式公式:其中,

(i=0,1,2,…，p)

若令t=xn+sh，則

(i=0,1,2,…，p)

當k，j，p取不同值時，得到不同類型的具體公式。如當k=1，j=0，p=0，1，2，…時，可得到阿達姆斯顯式

方法：

(i=0，1，2，…，p)

表8.4給出了顯式阿達姆斯公式的系數(shù)。在阿達姆斯顯式方法中，最常用的是p=3的情形：

又如當k=0，j=1，p=0，1，2，…時，得到阿達姆斯隱式公式(用n+1代替n)：

表8.5給出了隱式阿達姆斯公式的系數(shù)。在阿達姆斯隱式方法中，最常用的是p=3時的情形：

p=3時，阿達姆斯顯、隱式方法都是四階方法。

2.泰勒展開法(待定系數(shù)法)

用泰勒展開法構造線性多步法就是在線性多步法的一般公式

中，假設yn+j=y(xn+j)(j=0，1，2，…，k)，將y(xn+j)和y′(xn+j)在xn處用泰勒公式展開，與推導公式(8.13)式的方法完全類似，導出局部截斷誤差按h的升冪排列表達式，然后

再確定相應的待定系數(shù)αj，βj。

我們舉幾種特殊情形給予說明?？紤]線性兩步法，此時k=2，α2=1，則有

yn+2+α1yn+1+α0yn=h(β2fn+2+β1fn+1+β0fn)

由(8.14)式可得系數(shù)α0,α1,β0,β1,β2滿足:解此方程可得

從而一般線性二步方法可寫為

(8.15)對任意α0，有當α0≠－1時，c4≠0，(8.15)式是三階的；當α0=－1時，c4=0,c5≠0，(8.15)式就是二步四階公式(辛普森公式)：

其截斷誤差為又如當α0=－5時，方法是顯式的，否則就稱為隱式的，此時(8.15)式稱為二步四階米爾恩方法。

考慮線性三步法，此時k=3，α3=1，則有由(8.14)式可得系數(shù)α0,α1,α2,β0,β1，β2，β3滿足的方程為

此方程組中含有六個未知量、四個方程，因此有無窮多組解，也即有無窮多個三步方法。如令α0=α1=0，可得到

此時可得到阿達姆斯三階顯式方法:

其局部截斷誤差為又如,令α0=α1=0,α2=－1，且考慮隱式方法，可得到阿達姆斯四階隱式方法:

其截斷誤差為類似地，還可得到其他線性多步方法，如米爾恩方法:

哈明(Hamming)方法:

在應用線性多步法求解初值問題時，開始幾點處的函數(shù)值往往用別的方法計算。一般地選用與多步法同階的單步法，如龍格—庫塔方法、泰勒方法等。對線性隱式多步法，除開始幾點的函數(shù)值需單獨計算外，還需迭代求解或采用預估—校正方法求解。

例8.5應用四步四階阿達姆斯顯式方法求解初值問題：

取步長h=0.1。

解分析：應用四步顯式法必須有四個起步值，y0已知，而y1,y2，y3可用精度相同的四階龍格—庫塔方法求出。

步長h=0.1，節(jié)點xn=nh=0.1n(n=0,1,…，6)，由四階龍格—庫塔法算出：

y1=1.0048375,

y2=1.0187309,

y3=1.0408184

四步四階阿達姆斯顯式為已知fn=f(xn,yn)=xn－yn+1，h=0.1,xn=0.1n，則

由此算出

y4=1.0703231,

y5=1.1065356,

y6=1.1488186

例8.6以四步四階阿達姆斯顯式與隱式作為預估—校正公式解初值問題：

取步長h=0.1，并將結(jié)果與精確解比較。

解分析：由于預估公式與校正公式都是四階的，因此起步值也應按四階公式求出。

已知y0=1，按四階龍格—庫塔公式法算出：

y1=1.095446，y2=1.183217,

y3=1.264912

四階阿達姆斯顯式作為預估式(取步長h=0.1)，即四階阿達姆斯隱式作為校正式，即

該問題的精確解為。

此題的數(shù)值計算結(jié)果和精確結(jié)果及誤差列于表8.6。

8.5算法的穩(wěn)定性及收斂性

8.5.1算法的穩(wěn)定性

穩(wěn)定性在微分方程的數(shù)值解法中是一個非常重要的問題。因為微分方程初值問題的數(shù)值方法求解過程中存在著各種計算誤差，這些計算誤差，如舍入誤差等引起的擾動在傳播

過程中，可能會大量積累，對計算結(jié)果的準確性將產(chǎn)生影響，這就涉及到算法穩(wěn)定性問題。本節(jié)中僅介紹絕對穩(wěn)定性的概念。

定義8.5當在某節(jié)點xn上的yn值有大小為δ的擾動時，如果在其后的各節(jié)點xj(j>n)上的值yj產(chǎn)生的偏差都不大于δ，則稱這種方法是絕對穩(wěn)定的。

穩(wěn)定性不僅與算法有關，而且與方程中函數(shù)f(x,y)也有關，討論起來比較復雜。為簡單起見，通常只針對模型方程

y′=λy

(λ<0)

(8.16)

來討論。一般方程若局部線性化，也可化為上述形式。模型方程相對比較簡單，若一個數(shù)值方法對模型方程是穩(wěn)定的，并不能保證該方法對任何方程都穩(wěn)定，但若某方法對模型方程如此簡單的問題都不穩(wěn)定的話，也就很難用于其他方程的求解。下面以歐拉方法為例討論絕對穩(wěn)定性。

先考察向前歐拉方法的穩(wěn)定性。模型方程

y′=λy

(λ<0)

的歐拉公式為

yn+1=yn+λhyn=(1+λh)yn

設yn有誤差δn，則實際參與運算的，由此引起yn+1的誤差為δn+1，實際得到

即有

yn+1+δn+1=(1+λh)(yn+δn)從而有

δn+1=(1+λh)δn

要使|δn+1|<|δn|，必須有|1+λh|<1。因此向前歐拉方法的絕對穩(wěn)定域為|1+λh|<1。在復平面上，|1+λh|<1是以1為半徑，以－1為圓心的圓內(nèi)部，所以向前歐拉方法的絕對穩(wěn)

定域是圓域，如圖8-2所示。圖8-2向前歐拉方法的絕對穩(wěn)定域用向后歐拉方法解模型方程(8.16)式的計算公式為

yn+1=yn+λhyn+1

解出，則對于向后歐拉公式法，δn應滿足。由于λ<0，則有，故恒有|δn+1|<|δn|，因此向后歐拉方法是絕對穩(wěn)定的，它的絕對收斂域為|1－λh|>1。在復平面上，它是以1為中心，以1為半徑的圓外部，如圖8-3所示。圖8-3由圖8-2和圖8-3可知，向后歐拉方法的絕對穩(wěn)定域比向前歐拉方法的絕對穩(wěn)定域大得多，但可以證明這兩種方法的收斂階數(shù)是相同的，只是向前歐拉方法是顯式方法，向后歐拉方法是隱式方法。這也說明，隱式方法的穩(wěn)定性一般比同階的顯式方法的穩(wěn)定性要好得多。

用二級二階龍格—庫塔方法求解試驗方程(8.16)式的計算公式為

yn+1=yn+h［ω1λyn+ω2λ(yn+β21hyn)］

=［1+(ω1+ω2)λh+ω2β21(λh)2］yn利用二級二階龍格—庫塔方法的參數(shù)可知：

類似于向前歐拉方法的分析，可得二級二階龍格—庫塔方法的絕對穩(wěn)定域為

也就是|1+λh|<1，它與向前歐拉方法的絕對穩(wěn)定域相同。同理可得三級三階龍格—庫塔方法的絕對穩(wěn)定域是

絕對穩(wěn)定區(qū)間是(－2.51,0)。四級四階龍格—庫塔方法的絕對穩(wěn)定域是

絕對穩(wěn)定區(qū)間是(－2.78,0)。

二級二階、三級三階及四級四階龍格—庫塔方法的絕對穩(wěn)定域如圖8-4所示。圖8-48.5.2算法的收斂性

微分方程初值問題數(shù)值解法的基本思想是將微分方程轉(zhuǎn)化為差分方程來求解，并用計算值yn近似代替y(xn)。這種近似代替是否合理，還須看分割區(qū)間［xi－1，xi］的長度h越來越小時，即h=xi－xi－1→0時，yn→y(xn)是否成立。若成立，則稱該方法是收斂的；否則稱為不收斂。這里仍以歐拉方法為例，來分析收斂性。

歐拉格式如下：

yn+1=yn+hf(xn，yn)

設表示取yn=y(xn)時，按歐拉公式的計算結(jié)果，即

歐拉方法局部截斷誤差為設有常數(shù)，則

(8.17)

總體截斷誤差為

又

(8.18)由于f(x,y)關于y滿足李普希茲條件，即有

|f(xn,y(xn))－f(xn,yn)|≤L|y(xn)－yn|

(8.19)

將(8.19)式代入(8.18)式，有

再利用(8.17)式、(8.18)式可得到

(8.20)即

|εn|≤(1+hL)|εn－1|+ch2

上式反復遞推后，可得設xn－x0=nh≤T(T為常數(shù))，因為1+hL≤ehL，故(1+hL)n

≤enhL≤eTL把上式代入(8.20)式，得

若不計初值誤差，即ε0=0，則有

(8.21)

(8.21)式說明，當h→0時，εn→0，從而yn→y(xn)，所以歐拉方法是收斂的，且其收斂速度為O(h)，即具有一階收斂速度。

8.6一階常微分方程組與高階方程

前面已介紹了一階常微分方程初值問題的各種數(shù)值解法，對于一階常微分方程組，可類似得到各種解法，而高階常微分方程可轉(zhuǎn)化為一階常微分方程組來求解。8.6.1一階常微分方程組

對于一階常微分方程組的初值問題:

可以把單個方程y′=f(x,y)中的f,y看做向量來處理，這樣就可把前面介紹的各種算法推廣到求一階方程組初值問題中來。設xn=x0+nh(n=1,2,3,…)，yn，zn為節(jié)點xn上的近似解，則有下列求解格式：

(1)向前歐拉格式：

這是一步一階方法，為顯式。

(2)改進的歐拉格式：

這是一步二階方法。

(3)四階標準龍格—庫塔格式：

(8.22)其中，

(8.23)把節(jié)點xn上的yn和zn值代入(8.23)式，依次求出K1,L1,K2,L2,K3,L3,K4,L4，然后把它們代入(8.22)式，算出節(jié)點xn+1上的yn+1和zn+1值。

對于具有三個或三個以上方程的方程組的初值問題，也可用類似方法處理。此外，多步法也同樣可應用于求解方程組初值問題。

例8.7用改進的歐拉方法求解初值問題：

取步長h=0.1，保留6位小數(shù)。

解改進的歐拉公式為將f(xn,yn,zn)=xy－z,g(xn,yn,zn)=及h=0.1代入上式，得由初值y0=y(0)=1,z0=z(0)=2計算得8.6.2高階微分方程

高階微分方程(或方程組)的初值問題可以通過變量代換化為一階方程組來求解。例如，有二階微分方程的初值問題:

(8.24)作變換z=y′，則(8.24)式可化為一階方程組的初值問題:

(8.25)

對此問題就可用上小節(jié)的方法求解。此方法同樣可以用來處理三階或更高階的微分方程(或方程組)的初值問題。

例8.8采用改進的歐拉方法求解二階常微分方程初值問題:

(0≤x≤1)

取步長h=0.1，計算y(0.1)的近似值(最后結(jié)果保留小數(shù)點后

5位)。

解令z=y′，則所給微分方程初值問題等價于解該一階微分方程組，改進的歐拉公式為代入h=0.1,y0=－0.4,z0=－0.6,x0=0.0，計算得

由此可得

y(0.1)≈y1=－0.46200

例8.9應用四階龍格—庫塔方法求解二階常微分方程初值問題:

取步長h=0.1，計算y(0.1),z(0.1)的近似值(最后結(jié)果保留小數(shù)點后5位)。

解令z=y′，則所給微分方程初值問題等價于解該一階微分方程組，四階龍格—庫塔公式為其中,取步長h=0.1,x0=0,y0=0,z0=1,代入上式得則

由此可得

y(0.1)≈y1=0.10369,

z(0.1)≈z1=1.11043

8.7解微分方程的波形松弛方法

本節(jié)介紹波形松弛方法(也叫動力迭代方法)在求解微分方程的初值問題和邊值問題方面的應用。

通過積分方法可以求出一些特殊類型微分方程的解析解，但對維數(shù)較高而且形式復雜的微分方程，或者無法求得其解析解，或者所得到的解析解形式復雜難以得到所求解的性質(zhì)的問題，波形松弛方法是求解它們的一個有效工具。求解微分方程初值問題及邊值問題的迭代方法較多，本節(jié)僅介紹在畢卡(Picard)逐次逼近的基礎上發(fā)展起來的波形松弛方法在求解微分方程初值問題以及邊值問題方面的簡單應用。8.7.1微分方程初值問題的波形松弛方法

考慮微分方程的初值問題：

(8.26)

其中，f：R×Rn→Rn，y∈Rn,f滿足解的存在唯一性條件。為了敘述清楚，我們先以二階自治系統(tǒng)情形為例。

(8.27)

類似于線性代數(shù)方程組迭代解法中的迭代格式，有如下動力迭代格式。

1)雅可比動力迭代

雅可比動力迭代的基本迭代格式為

(k=0，1，2，…；x∈［x0,T］)

(8.28)通過逐次迭代，產(chǎn)生一系列函數(shù)序列{(y1(k)(x)，y2(k)(x))T},如果

則(y1(x)，y2(x))T就是(8.27)式的解。

雅可比動力迭代的優(yōu)點在于：對每一個k值，如果第k步已求出y1(k)(x)和y2(k)(x)，那么在進行第k+1步的求解過程中，(8.28)式中每個方程可以分別獨立地求解。這種分解格式特別對高維微分方程具有很大的優(yōu)勢，因為它可以進行并行計算，只在求出結(jié)果時，互相傳遞信息。

2)高斯—賽德爾動力迭代

高斯—賽德爾動力迭代的格式為

(k=0，1，2，…；x∈［x0,T］)

(8.29)同樣，經(jīng)過逐步迭代可得解序列{(y1(k)(x)，y2(k)(x))T}，當

時，就可得到(y1(x)，y2(x))T，它是(8.28)式的解。高斯—賽德爾動力迭代的優(yōu)點在于由(8.29)式中第一個方程求出解y1(k+1)(x)時，求解第二個方程時就用y1(k+1)(x)，而不用y1(k)(x)。這是因為在{(y1(k)(x),y2(k)(x))T}收斂時，y1(k+1)(x)比y1(k)(x)更接近于y1(x)，因此，在收斂情況下，我們把第一個方程已經(jīng)算出的y1(k+1)(x)代入第二個方程,同時也能使第一個方程的解序列收斂更快。高斯—賽德爾動力迭代的缺點在于只有當?shù)谝粋€方程求出y1(k+1)(x)時，才能計算第二個方程，因此不能在并行機上進行計算。例如,應用高斯—賽德爾動力迭代方法來計算下面兩個變量的微分方程初值問題的解：

初始函數(shù)選為，迭代一次有迭代兩次有

迭代三次有迭代n次有顯然，迭代序列{(y1(k)(t)，y2(k)(t))T}收斂到{(sint，cost)T},

即

y1(∞)(t)=sin(t)，y2(∞)(t)=cos(t)

事實上，我們?nèi)菀昨炞C該二階微分方程的解為

y1(t)=sin(t)，y2(t)=cos(t)

迭代序列收斂到精確解的過程如圖8-5和圖8-6所示，圖中給出了n=0，1，2，3，4，5時y1(n)(t)，y2(n)(t)的曲線圖。圖8-5圖8-6

3)SOR動力迭代

SOR動力迭代的基本格式為

(8.30)如果迭代函數(shù)序列收斂，即

則稱(y1(x)，y2(x))T為(8.27)式的解。

SOR迭代過程是把前一次得到的結(jié)果記為(z1(k)(x),z2(k)(x))T,然后把兩次所得結(jié)果進行加權平均，權系數(shù)為ω，所得結(jié)果作為下一次迭代的初值。SOR動力迭代的優(yōu)點是可以使迭代序列收斂速度更快，缺點是計算量大，同時權系數(shù)ω(也稱為松弛因子)的選取比較困難。對一般形式的n維非自治微分動力系統(tǒng)((8.26)式)，其動力迭代格式為

(8.31)此處，G和g是一般的分裂函數(shù)，且滿足：

G：［x0，T］×Rn×Rn×Rn→Rn，

G(x，y，y，y)=f(x，y)

g：［x0，T］×Rn×Rn→Rn，

g(x，y，y)=f(x，y)

事實上，(8.31)式的迭代格式包括了微分方程中畢卡逐次逼近的思想，它的一個更簡單的形式為

(8.32)此處，函數(shù)F滿足：

F：［x0，T］×Rn×Rn→Rn，

F(x，y，y)=f(x，y)

我們把(8.26)式中函數(shù)f的n個分量記為f1，f2，…，fn，且n維向量y(k)的n個分量記為y1(k)，y2(k)，…，yn(k)，則(8.26)式的雅可比、高斯—賽德爾、SOR動力迭代格式按分量寫出來分別為

(8.33)

(8.34)

(8.35)應當注意，在給出上述迭代格式以及迭代函數(shù)序列的收斂時沒有考慮收斂對自變量取值范圍的一致性問題。在應用動力迭代方法求解微分系統(tǒng)問題時一般要注意以下三個問題：

(1)積分區(qū)間上的誤差不是一致的。初始函數(shù)在起始階段能有比較好的逼近效果，但是隨著自變量x的增加，情況就會不斷惡化。

(2)每步迭代延長波形逼近真實解所需要的時間。

(3)計算波形和真實解之間的誤差未必在每個時間節(jié)點減少。從上一步迭代到下一步迭代，計算誤差反倒可能增加，尤其對于時間值x較大的范圍。8.7.2微分方程初值問題波形松弛方法的收斂問題

下面討論線性微分方程初值問題動力迭代的收斂性。收斂性結(jié)論主要基于壓縮映像原理。

線性微分方程初值問題為

(8.36)其中，A=(aij)n×n為已知矩陣；g(x)為n維向量函數(shù)，

g：R→Rn，y∈Rn,x∈R。對初值問題(8.36)式，其解析解為

而(8.36)式的動力迭代格式為

(8.37)其中，A=M－N。此處，M為對角矩陣或塊對角矩陣；N為非對角耦合矩陣。求解(8.37)式可得：

(8.38)

令

定義算子K為則(8.38)式可寫成不動點迭代形式：

y(k+1)(x)=Ky(k)(x)+Φ(x)(8.39)

容易證明(8.39)式的不動點是(8.36)式的解，且(8.36)式的解是(8.39)式的不動點。以y(x)表示(8.36)式的解，令ε(k)(x)=y(k)(x)－y(x)，則由(8.39)式可得

ε(k+1)(x)=Kε(k)(x)

從而有

ε(k)(x)=Kkε(0)(x)

(8.40)

此處Kk表示算子K的k重卷積。顯然,(8.40)式收斂的充分必要條件為ρ(K)<1。為了獲得第k步誤差ε(k)(x)對初始誤差ε(0)(x)的真實估計，首先給出如下定義。

定義8.6令‖·‖表示Cn中的任一種確定范數(shù)，則在［0,T］上，連續(xù)函數(shù)的最大范數(shù)定義為由定義8.6可以證明：

(8.41)

此處,表示算子γ的k重卷積γ*γ*…*γ，且γ*u(x)=

。由于已給T<+∞，因此有

‖γ‖T≤c

(8.42)且c=eT‖M‖‖N‖，則

由歸納可得

結(jié)合(8.39)式就有定理8.3。

定理8.3由動力迭代算法((8.37)式)以及K和γ(x)的定義，并由(8.42)式可得算子K的一個界為

(8.43)

注：由于，因此定理8.3表明對所有有限T值，由(8.37)式或(8.39)式產(chǎn)生的函數(shù)序列{y(k)(x)}收斂，收斂的速度可由K的分解表達式進行分析。由于

因此對任何固定T<+∞和非零λ，Rλ(K)是有界的。8.7.3微分方程邊值問題的波形松弛方法

微分方程邊值問題的基本形式為

(8.44)

對邊值問題的動力迭代，其基本格式完全類似于對初值問題的迭代。只要把迭代格式(8.33)式、(8.34)式和(8.35)式中相應的初值條件y(k+1)(x0)=y0改為y(k+1)(0)=y(k+1)(T)即可得到微分方程邊值問題的雅可比、高斯—賽德爾、SOR動力迭代格式。我們在此僅介紹微分方程邊值問題的動力迭代方法中的一個收斂性條件。

考慮微分方程邊值問題(8.44)式的動力迭代格式：

(x∈［0，T］；k=0，1，2，…)

(8.45)

其中，F(xiàn)：［0，T］×Rn×Rn→Rn，滿足F(x,y,y)=f(x,y)，

且x(0)(t)是一個初始猜測函數(shù)，滿足x(0)(0)=x(0)(T)。令‖·‖表示Rn中的2—范數(shù)，〈·，·〉表示標準內(nèi)積，即〈w，w〉=‖w‖2。假設F滿足下列條件：

(1)存在非負函數(shù)L1(x)，使得

‖F(xiàn)(x,u1,v)－F(x,u2,v)‖≤L1(x)‖u1－u2‖

(8.46)

其中，x∈［0,T］,v∈Rn,u1,u2∈Rn

。

(2)存在函數(shù)a(x)，使得

〈F(x,u,v1),F(x,u,v2)〉≤a(x)‖v1－v2‖

(8.47)

其中,x∈［0,T］,u∈Rn,v1,v2∈Rn,且。

先給出引理8.1。

引理8.1設a(x)∈R,Q(t)∈R,x∈［0,T］,如果φ(x)

滿足:

則有

證明首先證明當Q(x)≡0,x∈［0,T］時，φ0=0結(jié)論成立，因為

所以有即φ(x)≤0，x∈［0，T］。其次我們來證明其一般情形。

令

(x∈［0，T］)

則φ(x)滿足下列初值問題:我們定義

則有即

由前述情形可知，y(x)≤0，x∈［0，T］，即φ(x)≤ψ(x)，x∈［0，T］，從而完成了證明。

定理8.4對(8.44)式和(8.45)式，有

(8.48)

證明首先，令ε(k+1)(x)=y(k+1)(x)－y(x)，由(8.46)式和(8.47)式，有又由于

則對‖ε(k+1)(x)‖≠0,我們有

(8.49)

由引理8.1可知：

(8.50)此不等式對‖ε(k+1)(x)‖=0也成立。由于‖ε(k+1)‖(0)=

‖ε(k+1)‖(T)，因此由(8.50)式有

(8.51)

把(8.51)式代入(8.50)式，則有(8.48)式成立。

由引理8.1和定理8.4，可以得到定理8.5。

定理8.5對微分方程邊值問題(8.44)式和(8.45)式，如果函數(shù)f滿足條件(8.46)式和(8.47)式，且

則動力迭代序列{y(k)(x)}收斂到(8.44)式的解。

證明由于

定義范數(shù)‖·‖e如下：此處，‖·‖是Rn

中的范數(shù)，且λ是正常數(shù)，而

是有界的，因為此表達式內(nèi)的函數(shù)在［0，T］上是連續(xù)的，且

(8.52)

此處，因此我們有

(8.53)

此處

根據(jù)條件(8.52)式及定理假設條件，可以選取充分大λ，使L<1，即就是迭代格式(8.53)式收斂，從而{ε(k+1)(x)}收斂到零，證畢。下面我們就L1(x)=L1為常數(shù)，且a(x)=L2<0，且為常數(shù)給出一個數(shù)值計算的例子。

例8.10考慮如下二維系統(tǒng)：

此處，,我們?nèi)∪缦聞恿Φ袷?

(x∈［0,2π］；k=0,1,…)

表8.7給出了迭代次數(shù)與誤差的數(shù)值。8.8微分方程邊值問題的數(shù)值方法

由于高階微分方程經(jīng)過變量代換可以化為一階微分方程組，因此微分方程組邊值問題的一般形式為

其中，F(xiàn)：R×Rn→Rn，y∈Rn。假設F(x，y)滿足解的存在唯一性條件。在實際工程技術中經(jīng)常遇到二階微分方程：

y″=f(x,y,y′)

(x∈［a,b］)

(8.54)

為了確定唯一解，需要附加兩個定解條件。當定解條件取為解在區(qū)間［a,b］兩端點的狀態(tài)時，相應的定解條件稱為兩點邊值問題。

邊值條件有以下三類提法。

第一類邊界條件：

y(a)=α,

y(b)=β

當α=0或β=0時，稱為齊次邊界條件；否則稱為非齊次邊界條件。第二類邊界條件：

y′(a)=α,

y′(b)=β

當α=0或β=0時，稱為齊次邊界條件；否則稱為非齊次邊界條件。

第三類邊界條件：

y(a)－α0y′(a)=α1,

y(b)－β0y′(b)=β1

其中，α0≥0，β0≥0，α0+β0>0。當α1=0或β1=0時，稱為齊次邊界條件；否則稱為非齊次邊界條件。微分方程(8.54)式附加第一、二、三類邊界條件后，分別稱為第一、二、三類邊值問題。

下面僅介紹二階微分方程兩點邊值問題的打靶法和有限差分方法。8.8.1打靶方法

以二階微分方程第一類邊值問題為例來討論打靶法。它的基本原理是將邊值問題轉(zhuǎn)化為相應的初值問題：

(8.55)

來求解。令y1=y′，可以把(8.55)式化為

(8.56)至此，問題轉(zhuǎn)化為求適當?shù)膠，使初值問題(8.56)式的解y(x,z)在x=b處的值滿足右端邊界條件

y(b,z)=β

這樣,初值問題(8.56)式的解y(b,z)就是相應邊值問題的解。而求初值問題(8.56)式可用前面介紹的數(shù)值方法來求解，比如牛頓方法或其他迭代方法。先用線性插值方法在z=z0,

即y′(a)=z0

時求解初值問題(8.55)式，得解y(b,z0)=β0。若|β－β0|≤ε(ε為允許誤差)，則y(xj，z0)(j=0，1，2，…，m)是初值問題(8.56)式的數(shù)值解，也就是相應邊值問題的解。若|β－β0|>ε，則調(diào)整初始條件為y′(a)=z1，重新解初值問題(8.55)式，得解y(b,z1)=β1。若|β－β1|≤ε(ε為允許誤差)，則y(xj,z1)(j=0,1,2,…，m)為所求，否則再修正z1為z2