計算方法引論-第十一章

計算方法引論：

微分方程數(shù)值解法常微分方程初值問題的數(shù)值解法雙曲型方程的差分解法拋物型方程的差分解法橢圓型方程的差分解法有限元方法計算方法引論(第三版)第十一章常微分方程初值問題數(shù)值解法幾種簡單的數(shù)值解法R-K方法線性多步法預估—校正公式方程組和高階方程的數(shù)值解法自動選取步長的需要和事后估計剛性方程計算方法引論(第三版)引言問題求(11.1)的解y(x)在存在區(qū)間[a,b]中點列xn=xn–1+hn

(n=1,…,N)上的近似值yn. (11.1)hn＞0,步長.為敘述方便，假設(shè)hn不變,記為h.假設(shè)解y(x)在區(qū)間[a,b]上存在、唯一并充分光滑.

f(x,y)也充分光滑.

解法一步法：計算yn+1時只用到前一步的值.

多步法：計算yn+1要用到前k步的值.計算方法引論(第三版)

Euler方法

計算公式

(11.2)

由來幾何解釋過(xn,yn)沿斜率f(xn,yn)方向行一步

Taylor展開

數(shù)值微分

數(shù)值積分左矩形公式用于計算方法引論(第三版)

Euler方法算例

計算問題

y′=y-2x/y,0≤x≤1,y(0)=1，h=0.1計算結(jié)果

計算方法引論(第三版)Euler方法的誤差估計局部截斷誤差準確解代入

(11.2)式的剩余總體截斷誤差

,≤M2h2/2

計算方法引論(第三版)Euler方法的絕對穩(wěn)定性絕對穩(wěn)定性概念

用一個數(shù)值方法求解微分方程y′=λ

y,其中,λ是一個復常數(shù),對給定步長h>0,在計算yn時引入了誤差若這個誤差在后面的計算中所引起的誤差絕對值均不增加,就說這個數(shù)值方法對于這個步長h和復數(shù),λ是絕對穩(wěn)定的.Euler方法絕對穩(wěn)定性誤差方程:絕對穩(wěn)定條件:絕對穩(wěn)定區(qū)域:復平面上以λh=

-1為中心，1為半徑的圓

計算方法引論(第三版)

向后Euler方法

計算公式(11.10)

由數(shù)值積分右矩形公式導出誤差局部截斷誤差總體截斷誤差O(h)隱式方法迭代求解向后Euler方法是A穩(wěn)定的誤差方程:A穩(wěn)定:絕對穩(wěn)定區(qū)域為整個左半平面,Re(λh)<0

計算方法引論(第三版)

梯形方法

計算公式(11.14)

由數(shù)值積分梯形求積公式導出誤差局部截斷誤差總體截斷誤差O(h2)隱式方法迭代求解梯形方法是A穩(wěn)定的誤差方程:計算方法引論(第三版)?改進Euler方法

計算公式(11.17)

可視為一種預估—校正方法先用Euler方法預估pn+1=yn+hf(xn,yn)再用梯形方法校正yn+1=yn+h/2(f(xn,yn)+f(xn+1,pn+1))局部截斷誤差O(h3)總體截斷誤差O(h2)p階方法

總體截斷誤差為O(hp)的方法稱p階方法方法的總體截斷誤差為O(hp)則局部截斷誤差為O(hp+1)計算方法引論(第三版)?改進Euler方法算例

計算問題

y′=y-2x/y,0≤x≤1,y(0)=1，h=0.1計算結(jié)果

計算方法引論(第三版)

Taylor展開法一階方法到二次項得Euler方法yn+1=yn+hf(xn,yn)

r階方法到r+1次項得y′(x)=f(x,y(x))計算量很大.可作幾個點上f值的線性組合來構(gòu)造近似公式使其有r+1階局部截斷誤差,乃得R-K方法計算方法引論(第三版)幾個顯式R-K方法二階方法三階方法四階(經(jīng)典R-K)方法絕對穩(wěn)定區(qū)域誤差方程絕對穩(wěn)定區(qū)域

計算方法引論(第三版)

經(jīng)典R-K方法算例

計算問題

y′=y-2x/y,0≤x≤1,y(0)=1，h=0.2計算結(jié)果

計算方法引論(第三版)

Adams外插公式計算公式

(11.28)不能“自開始”可用四階R-K方法計算y1，y2，y3.由來出發(fā)點令F(x)=f(x,y(x))以xn,xn-1,xn-2,xn-3為節(jié)點作插值求積公式局部截斷誤差

計算方法引論(第三版)

Adams內(nèi)插公式計算公式

(11.30)不能“自開始”.

是隱式方法要迭代求解..由來出發(fā)點令F(x)=f(x,y(x))以xn+1,

xn,xn-1,xn為節(jié)點作插值求積公式局部截斷誤差

計算方法引論(第三版)Adams方法的穩(wěn)定性Adams外插法的絕對穩(wěn)定性誤差方程:通解：特征方程：絕對穩(wěn)定條件:絕對穩(wěn)定區(qū)域計算方法引論(第三版)Adams方法的穩(wěn)定性(續(xù))Adams內(nèi)插法的絕對穩(wěn)定性誤差方程:通解：特征方程：絕對穩(wěn)定條件:絕對穩(wěn)定區(qū)域計算方法引論(第三版)多步法穩(wěn)定性的注記多步法的絕對穩(wěn)定性本章從某步誤差對后面計算的影響考慮提出絕對穩(wěn)定性概念并討論了各方法的絕對穩(wěn)定性及絕對穩(wěn)定區(qū)域.一般多步法特征方程的根結(jié)構(gòu)復雜,絕對穩(wěn)定性也可定義為特征方程的根按模小于1(ρn→0,n→∞).穩(wěn)定性與收斂性尙有一些穩(wěn)定性概念以處理不同問題.類似于十二、三章有一種穩(wěn)定性保證(Dahlquist定理):

相容的差分方法(局部截斷誤差為:o(h))收斂的充要條件是其穩(wěn)定,即特征方程h=0時的根在單位圓,圓周上只有單根.因此本章各方法收斂計算方法引論(第三版)預估-校正公式顯式方法計算初值隱式方法進行迭代

Ｅuler-梯形Adams只迭代一次預估

校正

預估

校正

計算方法引論(第三版)?算例:R-K與AdamsP-C法計算問題

y′=y-2x/y,0≤x≤1,y(0)=1，h=0.1計算結(jié)果

計算方法引論(第三版)誤差估計Adams方法局部截斷誤差的估計由于得到估計據(jù)此可改善結(jié)果,乃有Adams預估—校正法的修正方案.

計算方法引論(第三版)修正的Adams預估-校正方法Adams預估—校正方法的一種修正方案原誤差的主要部分得以補償,是五階方法.y1，y2，y3另用單步法計算.并令

p3=c3=0.計算方法引論(第三版)常微分方程組解法常微分方程組n=2Euler方法向后Euler方法梯形公式k=0,1,…Adams外插公式計算方法引論(第三版)常微分方程組解法(續(xù))常微分方程組n=2Adams預估—校正方法k=0,1,…R-K方法計算方法引論(第三版)高階常微分方程解法化成方程組的形式求解二階方程化為一階方程組函數(shù)的特定形式計算公式可以化簡例如R-K公式

計算方法引論(第三版)事后估計誤差事后估計誤差依據(jù)是誤差的先驗估計方法一:利用同一公式不同步長的誤差來估計計算結(jié)果的誤差.數(shù)值積分中逐次分半法和Romberg求積用過.方法二:利用同階的不同公式的誤差來估計計算結(jié)果的誤差.本章四階Adams外插和內(nèi)插公式存在不同階的二公式,計算結(jié)果的差可估計其中一個的誤差.改善結(jié)果事后誤差估計可改善結(jié)果.Romberg求積公式Adams預估—校正法的修正方案事后誤差估計可調(diào)整步長.以一步r階公式為例介紹一種方法.計算方法引論(第三版)自動選取步長事后估計誤差局部截斷誤差為cnhr+1從xn–1出發(fā)兩步(步長為h)計算得從xn–1出發(fā)一步(步長為2h)計算得誤差可估計誤差并改善結(jié)果調(diào)整步長一種策略:如果

成立步長滿足要求.如過小,只1/2r容許誤差,步長加倍.否則步長減小一半自動變步長的R-K方法計算例子,效果顯著.尚有其它策略.多步法也可設(shè)計策略調(diào)整步長.不過調(diào)整代價較大.計算方法引論(第三版)

自動選取步長R-K方法算例

計算問題

其中

計算方法引論(第三版)

自動選取步長R-K方法算例

計算結(jié)果計算方法引論(第三版)

自動選取步長R-K方法算例

說明問題的解是以T=6.192169為周期的軌跡開始計算時給出步長h=0.4程序自動調(diào)節(jié)為0.0735892最小步長是第79步的0.000106687同樣精度要求用固定步長計算大約要計算62000步計算方法引論(第三版)Stiff問題Stiff現(xiàn)象

方程通解Reλ<0,|Reλ|愈大解衰減越快方程組二解衰減快慢差千倍系統(tǒng)解衰減速度局部地決定于Jacobi矩陣的特征值大小方程組求穩(wěn)態(tài)解(暫態(tài)解能忽略后)困難截斷誤差由慢變分量確定(方程組中y,對應(yīng)小特征值).穩(wěn)定性的要求決定于快變化分量(方程組中z,對應(yīng)大特征值).必須依穩(wěn)定性要求采用小步長,步數(shù)多,過程長.求穩(wěn)態(tài)解,