數(shù)值積分（基于MATLAB）課件 chapter5 常微分方程數(shù)值解；chapter6 非線性方程數(shù)值解

科學研究和工程實踐中，有很多實際問題的數(shù)學模型都是微分方程。利用微分方程理論，我們可以研究它們的一些性質(zhì)，對實際問題進行分析。但是，只有極少數(shù)特殊的方程有解析解。對于絕大部分的微分方程是沒有辦法求出它的解析解的。

常微分方程作為微分方程的基本類型之一，在自然界與工程界有很廣泛的應(yīng)用。很多問題的數(shù)學模型都可以歸結(jié)為常微分方程的定解問題。很多偏微分方程問題，也可以化為常微分方程問題來近似求解。常微分方程數(shù)值解本章討論常微分方程的數(shù)值解法考慮一階常微分方程的初值問題只要f(x,y)在[a,b]

R1上連續(xù)，且關(guān)于y

滿足Lipschitz

條件，即存在與x,y無關(guān)的常數(shù)L

使對任意x

[a,b]，和y1，y2

R1都有|f(x,y1)

f(x,y2)|

L|y1

y2|成立，則上述問題存在唯一解。

雖然有些微分方程初值問題的解是存在而且唯一的，但是，要求出其精確解y＝y(tǒng)(x)往往是非常困難的。為此，我們想法計算出解函數(shù)y(x)在一系列離散節(jié)點a=x0<x1<…<xn=b處的近似值yi

y

(xi)(i=1,2,…,n).——這就是微分方程數(shù)值解。節(jié)點間距hi=xi+1

xi為步長，通常采用等距節(jié)點，即取hi=h

(常數(shù))。求微分方程數(shù)值解的基本思想是：在這些節(jié)點上采用離散化方法，（通常用數(shù)值積分、微分、泰勒展開等）將初值問題化成關(guān)于離散變量的相應(yīng)問題。把這個相應(yīng)問題的解yn作為y(xn)的近似值。這樣求得的yn就是上述初值問題在節(jié)點xn上的數(shù)值解。一般說來，不同的離散化導致不同的方法。

微分方程數(shù)值解：計算微分方程解函數(shù)y(x)在節(jié)點a=x0<x1<…<xn=b處的近似值yi

y

(xi)(i=1,2,…,n).1.歐拉法：x0xi向前差商近似導數(shù)記為亦稱為歐拉折線法

x1定義在假設(shè)yi=y(xi)，即第

i

步計算是精確的前提下，考慮的截斷誤差Ri=y(xi+1)

yi+1稱為局部截斷誤差定義若某算法的局部截斷誤差為O(hp+1)，則稱該算法有p

階精度。微分方程數(shù)值解的精度

歐拉法的局部截斷誤差：=O(h2)=O(h1+1)歐拉法的精度為1階。例1

用歐拉法求初值問題在區(qū)間[0,0.10]上的數(shù)值解:當h=0.02時在區(qū)間[0,0.10]上的數(shù)值解。方程真解：nxnyny(xn)

n=y(xn)yn001.00001.0000010.020.98200.98250.000520.040.96500.96600.000530.060.94890.95030.001440.080.93360.93540.001850.100.91920.9230.0021當h=0.02時在區(qū)間[0,0.10]上的數(shù)值解。一階方程的初值問題與積分方程是等價的，當x=x1時，

借助于數(shù)值積分，求y(x1)的值

若用矩形公式2．梯形公式用梯形公式

則有——Euler公式于是有遞推格式：稱之為梯形公式。這是一個隱式的計算公式，欲求的yn+1需需要解一個方程。Euler公式是求微分方程數(shù)值解的很好的方法，它算法簡單，易于計算，但Euler方法有一個弱點就是誤差較大不能保證精度要求。梯形公式的弱點在于需要解一個方程。因此常采用將簡單的Euler公式與梯形公式相結(jié)合的方法：即預(yù)測-校正法：3．改進歐拉法——改進的Euler公式改進的Euler公式還有幾種不同形式：取步長h=0.1，分別用Euler方法和改進的Euler方法，求微分方程初值問題在區(qū)間[0，1.5]上的數(shù)值解。解：本題的精確解為，可用來檢驗近似解的精確程度。例題Euler方法：y0=1,yn+1=yn+hf

(xn,yn).y1=y(0.1)=y0+hf(x0,y0)=1+0.1

(120/1)=1.1,y2=y(0.2)=1.191818,……y11=y(1.0)=1.784778改進的Euler方法：y0=1,y1=y0+hf(x0,y0)=1.1,y1=1+01./2[(1-20/1)+(1.1-20.1/1.1]=1.095909,……y11=……y11=1.737869.xn歐拉法yn改進歐拉法yn準確解01110.11.11.0959091.0954450.21.1918181.1840961.1832160.31.2774381.2602011.2649110.41.3582131.3433601.3416410.51.4351331.4161021.4142140.61.5089661.4829561.4832400.71.5803381.5525151.5491930.81.6497831.6164761.6124520.91.7177791.6781681.6733201.01.7847701.7378691.7320511.11.851181.7958221.7888541.21.9174641.8522421.8439091.31.9840461.9073231.8973671.42.0514041.9612531.9493591.52.1200522.0142072.0000003．局部截斷誤差二元泰勒公式：設(shè)z=f(x,y)在點(x0，y0)的某一鄰域內(nèi)連續(xù)且直到有n+1階連續(xù)偏導數(shù)，(x0＋h，y0＋k)為此鄰域內(nèi)任一點，則有：歐拉法的截斷誤差：

改進歐拉法的截斷誤差：

建立高精度的單步遞推格式。單步遞推法的基本思想是從(xi,yi)點出發(fā)，以某一斜率沿直線達到(xi+1

,yi+1

)點。歐拉法及其各種變形所能達到的最高精度為2階?？疾旄倪M的歐拉法，可以將其改寫為：二、龍格-庫塔法首先希望能確定系數(shù)

1、

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

Step1:將K2在(xi,yi)

點作Taylor展開將改進歐拉法推廣為：),(),(][12122111phKyphxfKyxfKKKhyyiiiiii++==++=+llStep2:將K2代入第1式，得到Step3:將yi+1與y(xi+1)在xi點的泰勒展開作比較要求，則必須有：這里有個未知數(shù)，個方程。32存在無窮多個解。所有滿足上式的格式統(tǒng)稱為2階龍格-庫塔格式。取注意到此式就是改進的歐拉法。則有類似的有三階Runge-Kutta公式：問題:

為獲得更高的精度，應(yīng)該如何進一步推廣？其中

i

(i=1,…,m)，

i

(i=2,…,m)

和

ij

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

)

均為待定系數(shù)，確定這些系數(shù)的步驟與前面相似。)...,(......),(),(),(]...[1122112321313312122122111--++++++=+++=++==++++=mmmmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyybbbabbaballl

最常用為四級4階經(jīng)典龍格-庫塔法4階龍格――庫塔法截斷誤差階為O(h5)。取步長h=0.1，

用四階龍格――庫塔法解初值問題y

=x2

–y(0≤x≤1)y(0)=1的數(shù)值解

解：初值問題的精確解為y=x2-2x+2-e-

x例題用四階龍格――庫塔公式

xi

yi

y(xi)(精確解)0.000000,y0=1.000000y(x0)=1.0000000.100000,y1=0.905163y(x1)=0.9051630.200000,y2=0.821269y(x2)=0.8212690.300000,y3=0.749182y(x3)=0.7491820.400000,y4=0.689680y(x4)=0.6896800.500000,y5=0.643470y(x5)=0.6434690.600000,y6=0.611189y(x6)=0.6111880.700000,y7=0.593415y(x7)=0.5934150.800000,y8=0.590672y(x8)=0.5906710.900000,y9=0.603431y(x9)=0.6034301.000000,y10=0.632122y(x10)=0.632121寫成向量的形式：方程組和高階方程的數(shù)值解法各種方法都可以直接運用過來。以兩個方程的方程組為例Euler公式Runge-Kutta公式向量形式Runge-Kutta公式一般形式：例題：求微分方程組滿足初始條件在區(qū)間[0,1]上的數(shù)值解,取h=0.1.容易計算該方程組的精確解為：y=sinx

cosx+ex

z=sinx+cosx+ex

1x0=0.000000,y0=0.000000z0=1.000000x1=0.100000,y1=0.205342,yy1=0.210000,z1=1.210175,zz1=1.200009x2=0.200000,y2=0.422806,yy2=0.440006,z2=1.441955,zz2=1.400139x3=0.300000,y3=0.654718,yy3=0.690043,z3=1.697616,zz3=1.600716x4=0.400000,y4=0.903649,yy4=0.960182,z4=1.979671,zz4=1.802304x5=0.500000,y5=1.172443,yy5=1.250564,z5=2.290899,zz5=2.005729x6=0.600000,y6=1.464238,yy6=1.561426,z6=2.634373,zz6=2.212097x7=0.700000,y7=1.782506,yy7=1.893128,z7=3.013488,zz7=2.422813x8=0.800000,y8=2.131082,yy8=2.246191,z8=3.431996,zz8=2.639604x9=0.900000,y9=2.514206,yy9=2.621320,z9=3.894045,zz9=2.864540x10=1.000,y10=2.936564,yy10=3.019451,z10=4.404220,zz10=3.100055經(jīng)計算得數(shù)值解：1、2、確定方法，然后求解（0.202760.0881157）（0.2130070.0934037）（0.2237630.0988499）（0.2350520.104437）（0.2469020.110146）4階Runge-Kutta法，h=1高階方程則有：令用若干節(jié)點處的y

及y

值的線性組合來近似y(xn+1)。其通式可寫為：)...(...110111101knknnnknknnnffffhyyyy--+---+++++++++=bbbbaaa當

10時，為隱式公式;

1=0則為顯式公式。線性多步法（1）求出開頭幾個點上的近似值，即計算“表頭”；線性多步法：（2）利用逐步求后面點xk上的值yk。

基于數(shù)值積分的構(gòu)造法將在上積分，只要近似地算出右邊的積分，則可通過近似y(xn+1)。而選用不同近似式Ik，可得到不同的計算公式。得到若積分用節(jié)點作為積分點，則有積分系數(shù)這是顯格式，q+1階r+1步格式。r=max{p,q}局部截斷誤差為積分節(jié)點，可以構(gòu)造r+1步q+1階隱格式同樣，若以例：建立p=1,q=2的顯格式p=1，q=2，顯格式，積分區(qū)間為積分節(jié)點為所以例：建立p=2,q=2的隱格式p=2，q=2，隱格式，積分區(qū)間為積分節(jié)點為所以它的截斷誤差較顯格式小，通常也具有更好的穩(wěn)定性。1.阿當姆斯外推公式

以xn-2，xn-1，xn為節(jié)點作牛頓向后插值多項式P2(x)。其中對于一般的差分方程仍然考慮最簡單的模型，即只有初值產(chǎn)生誤差，看看這個誤差的傳播。差分方程的絕對穩(wěn)定性由初始誤差產(chǎn)生了差分解的誤差，實際上是同一差分方程，取不同初值所得到的2組差分解之間的差。這個差不僅于差分方程本身有關(guān)，而且與微分方程本身有關(guān)。如果微分方程本身是不穩(wěn)定，那就沒理由要求這2組解充分接近。因此，差分方程的穩(wěn)定性概念是建立在微分方程穩(wěn)定的基礎(chǔ)上的。把這個典型微分方程規(guī)定為：差分方程運用到如上的微分方程后，可以得到對于給定的初始誤差，誤差方程具有一樣的形式定義：差分方程稱為絕對穩(wěn)定的，若差分方程作用到微分方程時，對任意的初值，總存在左半復(fù)平面上的一個區(qū)域，當在這個區(qū)域時，差分方程的解趨于0。這個區(qū)域稱為穩(wěn)定區(qū)域數(shù)值解的穩(wěn)定性例：考察初值問題在區(qū)間[0,0.5]上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。0.00.10.20.30.40.5精確解改進歐拉法

歐拉隱式歐拉顯式

節(jié)點xi

1.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.05901071.0000001.3750001.8906252.5996103.5744644.914889龍格－庫塔例：向后Euler公式的穩(wěn)定性誤差方程：210ReImg考察隱式歐拉法可見絕對穩(wěn)定區(qū)域為：210ReImg注：一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式法的好。3階Runge－Kutta顯式

1~4階方法的絕對穩(wěn)定區(qū)域為k=1k=2k=3k=4-1-2-3---123ReImg插值公式的余項為則積分公式的截斷誤差為k=3時的外推公式為余項為：將差分表示成函數(shù)值的和的形式：二階阿當姆斯外推公式可改寫為：三階阿當姆斯外推公式可改寫為：2.阿當姆斯內(nèi)插公式將被積函數(shù)用以xn-1，xn，xn+1為插值節(jié)點的內(nèi)插多項式得到：k=1k=2時阿當姆斯外推法與內(nèi)插聯(lián)合起來考慮常微分方程的邊值問題：

其中p(x)，q(x)和f(x)均為[a,b]上給定的函數(shù)，

，

為已知數(shù)。假定p(x)、q(x)及f(x)均為[a,b]上充分光滑的函數(shù)，且q(x)≤0，這時，邊值問題存在連續(xù)可微的解，且唯一。解二階常微分方程邊值問題的差分法用差分法解邊值問題的主要步驟是：（1）將區(qū)間[a,b]離散化；（2）在這些節(jié)點上，將導數(shù)差商化，從而把微分方程化為差分方程；(3)解差分方程――實際上就是解線代數(shù)方程組。將[a,b]區(qū)間用節(jié)點分成N等分，其中x0=a與xN=b稱為邊界點，而x1,x2,…,xN-1稱為內(nèi)點。例9.7

試用差分法解方程解

將[0,1]劃分為四等分，即取，得五個節(jié)點差分方程為將它改寫成

在每個內(nèi)點列方程得由追趕法公式解得：y3=1.4855y2=1.2802y1=0.7753二分法迭代法迭代法的加速（Aitken加速法、Steffensen迭代法）牛頓迭代法第6章非線性方程的數(shù)值解§1.非線性方程實根的對分法(二分法)2）對[a1,b1]重復(fù)1)的計算，并產(chǎn)生[a2,b2]二分法的收斂性ax*x0ba1b1§2.迭代法迭代過程的幾何表示