數(shù)值分析(常微分方程數(shù)值解)

解一階常微分方程歐拉法局部截斷誤差與p階精度Range-Kutta公式一階常微分方程組和二階方程線性多步法簡介《數(shù)值分析》23*常微分方程(ODE)的定解問題主要有初值問題(IVP)和邊值問題(BVP)兩大類。

急性傳染病數(shù)學(xué)模型(SIR模型)常微分方程的初值問題是描述系統(tǒng)發(fā)展演變的重要工具。洛倫茲吸引子

lorenzgui關(guān)于愛情的動力系統(tǒng)

Loveraffairesanddifferentialequations*一階常微分方程初值問題:離散化(discretization)

取定離散點x0<x1<x2<···<xN,步長h

其中y=y(x)是未知函數(shù),右端函數(shù)f(x,y)是已知函數(shù),初值y0是已知數(shù)據(jù)。求未知函數(shù)y(x)在離散點處的近似值y1,y2,y3,·····,yN考察點xn列出微分方程Euler法與修正的Euler法其中yn是y(xn)的近似。足夠多的點來近似連續(xù)對象*求解常微分方程初值問題的Euler方法

取定步長:h,記

xn=x0+nh,(n=1,2,···,N)稱計算格式:yn+1=yn+hf(xn,yn)為顯式Euler公式。對應(yīng)的求初值問題數(shù)值解的方法稱為Euler方法。例1.

用Euler法求初值問題的數(shù)值解。解:記

f(x,y)=y－xy2,xn=nh(n=0,1,2,···,N)

由Euler公式得:yn+1=yn+h(yn－xnyn2)(n=0,1,···,N)*取步長

h=2/10,2/20,2/30,2/40,用Euler法求解的數(shù)值實驗結(jié)果如下.N10203040h0.20.10.06670.05誤差0.10590.05210.03420.0256解析解:o——數(shù)值解----——準(zhǔn)確解

*y'=f(x,y)左矩形公式用數(shù)值積分方法離散化常微分方程顯式Euler法求解常微分方程初值問題的數(shù)值方法可以由很多,各種方法都可以看成使用不同的數(shù)值積分方法計算*梯形公式:

右矩陣公式:

隱式Euler法*預(yù)-校方法算法如下k1=f(xn,yn),

k2=f(xn+1,yn+hk1)

,預(yù)報-校正方法(修正的Euler法):預(yù)報校正*預(yù)-校方法(h=0.2)誤差最大值0.0123n10203040h0.20.10.06670.05誤差20.01230.00260.00115.9612e-004誤差10.10590.05210.03420.0256歐拉方法(h=0.2)誤差最大值0.1059*

基本概念單步法:即為了求得后一點xn+1上數(shù)值解yn+1,只要知道一點xn的數(shù)值解yn就可以了。多步法:即計算yn+1時,要用到前面多點的數(shù)值解yn,yn-1,···,yn-r等信息。顯式格式:數(shù)值解yn+1可以用xn和yn

解析表出。這類格式稱為顯式格式。反之稱為隱式格式。對隱式格式而言,計算數(shù)值解yn+1形式上需要解方程(組)。顯式方法易于計算,但通常數(shù)值穩(wěn)定性較差。而隱式方法通常穩(wěn)定性較好,但計算上顯然很不方便。將顯式方法和隱式方法聯(lián)合使用是重要的思想方法。*

Tn+1=y(xn+1)-yn+1稱為局部截斷誤差。設(shè)

yn=y(xn)即由Taylor公式Euler公式:yn+1=yn+hf(xn,yn)的局部截斷誤差Euler公式的局部截斷誤差記為

O(h2),則稱Euler公式具有1階精度。局部截斷誤差僅考慮xn到xn+1的局部情況,并假定xn之前的計算沒有誤差*若局部截斷誤差為O(hp

+1)則稱顯式單步法具有

p階精度。例2.

證明梯形方法具有2階精度從x0開始計算,如果考慮每一步產(chǎn)生的誤差直到xn,則誤差y(xn)-yn稱為數(shù)值方法在節(jié)點xn處的整體誤差。粗略的講,整體誤差把許多步的局部截斷誤差加在一起,而步數(shù)正比于1/h(即步長倒數(shù)),因此整體誤差O(hp)。*

(Runge-Kutta)龍格-庫塔法改進(jìn)的Euler方法這樣理解:它用xn和xn+1兩個點的斜率值K1和K2取算術(shù)平均作為平均斜率,而xn+1處的斜率值則利用已知信息xn來預(yù)報。RK方法是一大類的方法,其基本思想是采用如下形式:是否可以推廣改進(jìn)的Euler方法？*取點個數(shù)s稱為龍格-庫塔法的級數(shù)。一般的RK公式如下給出:

RK方法是單步法(為了求得后一點xn+1上數(shù)值解yn+1,只要知道一點xn的數(shù)值解yn)。我記得我的朋友JohnvonNeumann曾經(jīng)說過,用四個參數(shù)我可以擬合出一頭大象,而用五個參數(shù)我可以讓它的鼻子擺動-費米下面以2級方法為例子具體介紹龍格-庫塔法。構(gòu)造的基本思想是選擇適當(dāng)?shù)南禂?shù)使得方法的局部截斷誤差階數(shù)盡可能高。**三階Range-Kutta公式一般形式y(tǒng)n+1=yn+h[k1+4k2+k3]/6k1=f(xn,yn),k2=f(xn+0.5h,yn+0.5hk1)k3=f(xn+h,yn–hk1+2hk2)四階Range-Kutta公式一般形式y(tǒng)n+1=yn+h[k1+2k2+2k3+k4]/6k1=f(xn,yn),k2=f(xn+0.5h,yn+0.5hk1)k3=f(xn+0.5h,yn+0.5hk2),k4=f(xn+h,yn+hk3)*例4數(shù)值實驗:幾種不同求數(shù)值解公式的誤差比較n10203040h0.20.10.06670.05RK46.862e-0053.747e-0067.071e-0072.186e-007RK30.00121.529e-0044.517e-0051.906e-005RK20.01230.00260.00115.9612e-004Euler0.10590.05210.03420.0256*常微分方程初值問題的MATLAB數(shù)值求解命令:[x,y]=ode23('f',[a,b],y0)f=inline('y-x.*y.^2');[x,y]=ode23(f,[0,2],1)符號求解命令:dsolve('eqn1',...)symsxydsolve('Dy=y-x*y^2','y(0)=1','x')解析解:參考:TheMatlabODEsuite,SIAMJScientificComputing*一階常微分方程組歐拉公式:(n=0,1,2,·······,N)一階常微分方程組向量形式記*修改的歐拉公式:(n=0,1,·······,N)經(jīng)典龍格-庫塔公式:*洛倫茲模型取=8/3，=10，=28。x(0)=0，y(0)=0，z(0)=0.01。t∈[0，80]，微分方程右端函數(shù):