劉禹希， 王 密

(武漢大學(xué)a.數(shù)學(xué)與統(tǒng)計學(xué)院;b.遙感信息工程學(xué)院，湖北 武漢 430061)

1 Cauchy問題的一種數(shù)值解法

Cauchy問題的求解是常微分方程理論中十分重要的課題，也是一階常微分方程初值問題中最基礎(chǔ)的一環(huán)。關(guān)于這種初值問題解的存在性理論，有以下著名的Picard定理：

定理1(Picard定理)

設(shè)初值問題：

(1)

式(1)中f(x,y)在矩形區(qū)域

R:|x-x0|≤a,|y-y0|≤b

(2)

內(nèi)連續(xù)，且對y滿足Lipshitz條件，則(E)在區(qū)間I=[x0-h,x0+h]上有且僅有一個連續(xù)可微解，其中常數(shù)

本文僅研究形如(E)，滿足條件和Lipshitz條件的常微分方程初值問題的數(shù)值解。

在該定理證明過程中將微分方程通過代數(shù)變形轉(zhuǎn)化為積分方程，并構(gòu)造了如下Picard序列：

(3)

式(3)中y(x0)=y0

事實上通過選定初值并對該序列進行迭代，經(jīng)過多次操作便可得到該微分方程的一個數(shù)值解，其截斷誤差有估計:

(4)

將以上數(shù)值求解方法稱為Picard方法。

定理 1表明：對一類特殊的Cauchy問題,其解一定可以唯一地定義在一個足夠小的閉區(qū)間中。但是對于整體積分曲線的研究，仍需要研究其局部解的最大存在區(qū)間問題，于是便有如下的微分方程解的延拓定理。

以下不加證明地給出該定理：

定理2(解的延拓定理)

設(shè)p為區(qū)域G內(nèi)一點，設(shè)T為微分方程經(jīng)過p的任一積分曲線。則積分曲線T將在區(qū)域G內(nèi)延伸到邊界。

2 數(shù)值積分方法

解決簡單數(shù)值積分的常見方法有：利用牛頓—柯斯特公式求解，利用復(fù)合求積法得到復(fù)合梯形公式和復(fù)合辛普森公式求解。前面介紹的求積方法可提高求積精度，如若精度仍不夠，則可通過將步長逐次減半的方式來提高精度,將改進后的公式稱之為龍貝格求積公式。

貫穿于數(shù)值積分中最經(jīng)典的思想便是“分割求和”，于是有以下著名的復(fù)合梯形公式：

公式1(復(fù)合梯形公式)

將[a,b]分為n等份，其節(jié)點為xi=a+ih,i=0,1,2,…,n,h=b-a/n，則有：

(5)

如對復(fù)合梯形公式進行代數(shù)變形即可導(dǎo)出其遞推公式(6):

(6)

式(6)中T2n表示在Tn基礎(chǔ)上步長h減半后的積分值。

為了便于刻畫復(fù)合梯形公式中的積分項，可以采用冪級數(shù)展開的方法逼近計算，于是便有著名的理查德森定理。

以下不加證明地給出該定理：

定理3(理查德森定理)

設(shè)f(x)∈C∞[a,b]，則有：

T(h)=I+α1h2+α2h4+…+αlh2l+…

(7)

其中T(h)=Tn,系數(shù)αl(l=1,2,…)與h無關(guān)。

從式(7)可以看出，用近似積分值I，誤差階為O(h4)，誤差階得到了提高。類似地，將計算I的近似值的誤差階由O(h2)提高到O(h4)的方法稱為外推算法，也稱為理查德森外推算法[3]。

以上是數(shù)值分析中一個重要的技巧，當(dāng)真值與近似值的誤差能夠表示成h的冪級數(shù)時，都可以使用外推算法來提高精度。以下將充分利用這種方法，進行適當(dāng)?shù)母倪M，推導(dǎo)出高精度的龍貝格求積算法。

3 龍貝格求積算法

龍貝格算法是在積分區(qū)間逐次分半的過程中，對用復(fù)化梯形法產(chǎn)生的近似值進行加權(quán)平均，以獲得準(zhǔn)確度較高的一種方法。具有公式簡練，使用方便，結(jié)果更可靠的特點。

事實上，對于I的計算可以轉(zhuǎn)化為對于近似迭代序列{Tn}的計算，通過迭代序列各項的比對，經(jīng)過代數(shù)變形整理即可。

龍貝格算法推導(dǎo)如下：

設(shè)m次外推加速為式(8)：

(8)

余項為式(9)：

Tm(h)=I+δ1h2(m+1)+δ2h2(m+2)+…

(9)

(10)

式(10)稱為龍貝格求積算法[4]，算法運行過程如下：

3.求加速值。

4 算法應(yīng)用實例

下面以一個此類型初值問題的數(shù)值求解來展示這種方法的實際應(yīng)用和重要意義。

4.1 例子

(11)

這是一類標(biāo)準(zhǔn)的一階線性的微分方程初值問題。由定理1，可知該初值問題存在唯一解；由定理2，可知該解可以延拓至整個區(qū)間[0,1]上成為一條積分曲線。

同時以上常系數(shù)線性微分方程可以通過求解其特征方程得到特解，并用定系數(shù)求得通解的方法得到其解析解為：y=ex+2x-x2.

4.2 Picard迭代算法

根據(jù)以上定理以及算法的鋪墊，整理可得常微分方程數(shù)值求解中的新型Picard算法，流程如下：

(1)寫出Picard迭代序列：

(3)將以上數(shù)值積分項代入Picard序列并計算出yn+1(x)；

(4)交替進行2和3直至迭代序列滿足|yn+1(x)-yn(x)|≤0.01，停止迭代。

4.3 結(jié)果展示

以下將通過比對常微分方程數(shù)值求解中著名的梯形迭代格式算法來體現(xiàn)Picard迭代算法的優(yōu)越性[5]。

利用MATLAB對以上Picard算法進行迭代，可得如下結(jié)果：

由圖1可知：該序列僅僅迭代四次便可很好地收斂到解析解，收斂速度非?？欤@體現(xiàn)了Picard迭代算法的高效性。

圖1 迭代序列

以下是誤差表格1。

表1 Picard迭代誤差表

將三種迭代格式的結(jié)果用MATLAB作圖呈現(xiàn)：

由圖2可知，相比于常微分方程數(shù)值求解中的梯形迭代格式，本文構(gòu)建的Picard迭代序列有著更好的收斂精度。

圖2 算法比較

分別取x=0.2,0.4,0.6,0.8,1.0，將Picard序列和梯形迭代序列與標(biāo)準(zhǔn)解析解的誤差做成誤差表如下表2。

表2 不同數(shù)值求解方法誤差表

由表2可知，當(dāng)變量x增大時，Picard迭代誤差及梯形迭代誤差均增加，但前者增加的速度明顯慢于后者。因此相比于傳統(tǒng)數(shù)值迭代方法，新型Picard迭代方法的效率明顯較高。

事實上，在面對一類特殊的常微分方程初值問題時，相比于傳統(tǒng)迭代方法，往往用Picard迭代方法可以更快地收斂到解，這也是迭代不動點方法的重要作用之一。

5 結(jié) 語

對傳統(tǒng)Picard迭代算法進行改進，引入復(fù)合龍貝格數(shù)值積分方法計算迭代序列，將求解精度進一步優(yōu)化，并結(jié)合實例驗證其迭代穩(wěn)定性同時和傳統(tǒng)數(shù)值方法進行對比，得到了很好的結(jié)果；同時，該算法也體現(xiàn)出迭代不動點數(shù)值方法極快的收斂速度，在面對更復(fù)雜的多項式型Cauchy問題時，該方法的求解所占內(nèi)存相對較小，有著很好的應(yīng)用前景。但是在面對一系列略病態(tài)的Cauchy問題時，該方法求龍貝格數(shù)值積分部分的運算復(fù)雜度仍然會很大，如何針對性地尋找更便捷的方法也是未來常微分方程數(shù)值求解方法應(yīng)當(dāng)著重研究的一個領(lǐng)域。