常微分方程的數(shù)值解法課件

第六章

常微分方程的數(shù)值解法.本章內(nèi)容§6.1引言§6.2歐拉方法§6.3龍格—庫塔方法§6.4邊值問題的數(shù)值方法.一.問題提出有一個或多個導數(shù)及其函數(shù)的方程式稱為微分方程，在工程中常遇到求解微分方程的問題?！?.1引言.§6.1引言二.兩類定解問題◆常微分方程的定解問題有兩種基本類型類:初值問題和邊值問題◆定解指已知因變量和/或其導數(shù)在某些點上是已知的(約束條件)?！?.邊值問題約束條件為已知，在自變量的任一非初值上，已知函數(shù)值和/或其導數(shù)值，如例如,受連續(xù)分布橫向荷載的變截面簡支梁彎曲問題.q(x)xwOl●2.初值問題例如，單自由度系統(tǒng)的非線性受迫振動.實際問題中還存在初邊值混合問題，如梁在橫向激勵下的彎曲振動。高階常微分方程可以化成一階的常微分方程組.很多微分方程的解不能用初等函數(shù)來表示，有時即使能夠用解析式表示其解，但計算量太大而不實用(表達式過于復雜)。需要用數(shù)值方法來求解，一般只要求得到若干個點上的近似值或者解的簡單的近似表達式(精度要求滿足即可)?！?.1引言.§6.1引言.§6.1引言.初值問題的常見解法單步法：利用前一個單步的信息(一個點)，在y=f(x)上找下一點yi，有歐拉法，龍格－庫格法。預測校正法：多步法，利用一個以上的前點信息求f(x)上的下一個yi，常用迭代法，如改進歐拉法，阿當姆斯法?！?.1引言.§6.2歐拉方法及其改進

Euler’sMethod

內(nèi)容一.歐拉格式二.Euler預估—校正法三.誤差估計、收斂性和穩(wěn)定性.6.2.1歐拉公式：/*Euler’sMethod*/向前差商近似導數(shù)記為§6.2歐拉方法及其改進.xP0P1P2P3P4PnyO.§6.2歐拉方法及其改進.§7.2歐拉方法hxiyi真值y(xi)誤差y(xi)-yih=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.00227.§6.2歐拉方法及其改進.§7.2歐拉方法§6.2歐拉方法及其改進.§7.2歐拉方法.§6.2歐拉方法及其改進6.2.3隱式歐拉法/*implicitEulermethod*/向后差商來近似導數(shù))1,...,0(),(111-=+=+++niyxfhyyiiii由于未知數(shù)yi+1

同時出現(xiàn)在等式的兩邊，不能直接得到，故稱為隱式

/*implicit*/

歐拉公式，而前者稱為顯式

/*explicit*/歐拉公式。一般先用顯式計算一個初值，再用隱式法（迭代）求解。.§6.2歐拉方法及其改進.§6.2歐拉方法及其改進.顯、隱式兩種算法的平均。需要迭代求解，能否不迭代？§6.2歐拉方法及其改進6.2.4梯形格式.

§6.2歐拉方法及其改進/*predictor-correctormethod*/Step1:先用顯式歐拉公式作預測，算出),(1iiiiyxfhyy+=+Step2:再將代入隱式梯形公式的右邊作校正，得到1+iy)],(),([2111+++++=iiiiiiyxfyxfhyy.§7.2歐拉方法及其改進.§7.2歐拉方法.§7.2歐拉方法ixiEuler方法yi改進Euler法yi精確解y(xi)01234567891000.10.20.30.40.50.60.70.80.9111.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.7177791.78477011.0959091.1840961.2662011.3433601.4164021.4859561.5525151.6164761.6781681.73786911.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051計算結果如下:.§6.3龍格—庫塔方法內(nèi)容一.2階龍格—庫塔格式三.高階龍格—庫塔格式.單步法：即利用前一個節(jié)點的函數(shù)值yi，計算后一個節(jié)點的函數(shù)值yi+1。目的：建立高精度的單步遞推格式。單步遞推法的基本思想是從(xi,yi)點出發(fā)，以某一斜率沿直線達到(xi+1,yi+1)點。歐拉法及其各種變形所能達到的最高精度為2階?！?.3龍格—庫塔方法二.2階龍格—庫塔格式.斜率一定取K1K2的平均值嗎？步長一定是一個h嗎？§6.3龍格—庫塔方法.§7.3龍格—庫塔方法首先希望能確定系數(shù)

1、

2、p，使得到的算法格式有2階精度，即在的前提假設下，使得

Step1:

將K2在(xi,yi)點作Taylor展開將改進歐拉法推廣為：),(),(][12122111phKyphxfKyxfKKKhyyiiiiii++==++=+ll.Step3:

將y(xi+1)在xi點的泰勒展開并與yi+1作比較要求，則必須有：這里有3個未知數(shù)，2個方程。存在無窮多個解。所有滿足上式的格式統(tǒng)稱為2階龍格-庫塔格式。注意到，就是改進的歐拉法。Step2:

將K2代入第1式，得到.其中

i

(i=1,…,m)，

i

(i=2,…,m)

和

ij

(i=2,…,m;j=1,…,i1

)

均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。)...,(......),(),(),(]...[1122112321313312122122111--++++++=+++=++==++++=mmmmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyybbbabbaballl問題:為獲得更高的精度，應該如何進一步推廣？§6.3龍格—庫塔方法三.高階龍格—庫塔格式.3階龍格-庫塔法§6.3龍格—庫塔方法.§6.3龍格—庫塔方法最常用為四級4階經(jīng)典龍格-庫塔法.注：

龍格-庫塔法的主要運算在于計算Ki

的值，即計算f

的值。Butcher于1965年給出了計算量與可達到的最高精度階數(shù)的關系：753可達到的最高精度642每步須算Ki的個數(shù)

由于龍格-庫塔法的導出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用低階算法而將步長h取小?！?.3龍格—庫塔方法.§7.3龍格—庫塔方法例：使用高階R-K方法計算初值問題解:（1）使用三階R-K方法(公式)§6.3龍格—庫塔方法.§7.3龍格—庫塔方法其余結果如下:

ixiK1K2K3

yi1.00000.10001.00001.10251.25551.11112.00000.20001.23451.37551.59451.24993.00000.30001.56241.76372.09221.42844.00000.40002.04042.34232.86581.66645.00000.50002.77683.25874.16341.9993§6.3龍格—庫塔方法.§7.3龍格—庫塔方法（2）如果使用四階R-K方法(公式)§6.3龍格—庫塔方法.§7.3龍格—庫塔方法其余結果如下:

ixiK1K2K3K4yi1.00000.10001.00001.10251.11331.23511.11112.00000.20001.23461.37561.39211.56331.25003.00000.30001.56251.76391.79082.04231.42864.00000.40002.04082.34282.38922.78051.66675.00000.50002.77773.26003.34764.00572.0000§6.3龍格—庫塔方法.步長過大，達不到精度要求；步長過小，雖然局部截斷誤差小，加大了計算工作量，舍入誤差的累積增大。解決途徑之一——引入變步長技術，常用的有Richardson外推法。從結點xi出發(fā)，先以h為步長，通過一步計算出y(xi+1)的近似值§6.4步長的自動選擇....§7.2歐拉方法/*Convergency*/§6.5收斂性與穩(wěn)定性.例：就初值問題考察歐拉顯式格式的收斂性。解：該問題的精確解為歐拉公式為對任意固定的x

=xi=ih，有

§6.5收斂性與穩(wěn)定性.§7.2歐拉方法例：考察初值問題在區(qū)間[0,0.5]上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。0.00.10.20.30.40.5精確解改進歐拉法

歐拉隱式歐拉顯式

節(jié)點xi1.0000

2.00004.0000

8.00001.6000101

3.2000101

1.00002.5000101

6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107原因??!取步長

h=0.1§6.5收斂性與穩(wěn)定性.§7.2歐拉方法§6.5收斂性與穩(wěn)定性.一般分析時為簡單起見，只考慮試驗方程常數(shù)，可以是復數(shù)§6.5收斂性與穩(wěn)定性.§6.5收斂性與穩(wěn)定性.Euler法的絕對穩(wěn)定區(qū)域.§7.3龍格—庫塔方法例：隱式龍格-庫塔法而顯式1~4階方法的絕對穩(wěn)定區(qū)域為其中2階方法的絕對穩(wěn)定區(qū)域為0ReImk=1k=2k=3k=4-1-2-3---123ReIm無條件穩(wěn)定.§6.6一階常微分方程組與高階方程

可以把單個方程中的