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

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

lorenzgui關于愛情的動力系統

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

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

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

取定步長:h,記

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

用Euler法求初值問題的數值解。解:記

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法求解的數值實驗結果如下.N10203040h0.20.10.06670.05誤差0.10590.05210.03420.0256解析解:o——數值解----——準確解

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

右矩陣公式:

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

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

,預報-校正方法(修正的Euler法):預報校正*預-校方法(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上數值解yn+1,只要知道一點xn的數值解yn就可以了。多步法:即計算yn+1時,要用到前面多點的數值解yn,yn-1,···,yn-r等信息。顯式格式:數值解yn+1可以用xn和yn

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

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

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開始計算,如果考慮每一步產生的誤差直到xn,則誤差y(xn)-yn稱為數值方法在節點xn處的整體誤差。粗略的講,整體誤差把許多步的局部截斷誤差加在一起,而步數正比于1/h(即步長倒數),因此整體誤差O(hp)。*

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

RK方法是單步法(為了求得后一點xn+1上數值解yn+1,只要知道一點xn的數值解yn)。我記得我的朋友JohnvonNeumann曾經說過,用四個參數我可以擬合出一頭大象,而用五個參數我可以讓它的鼻子擺動-費米下面以2級方法為例子具體介紹龍格-庫塔法。構造的基本思想是選擇適當的系數使得方法的局部截斷誤差階數盡可能高。**三階Range-Kutta公式一般形式yn+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公式一般形式yn+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數值實驗:幾種不同求數值解公式的誤差比較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數值求解命令:[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)經典龍格-庫塔公式:*洛倫茲模型取=8/3，=10，=28。x(0)=0，y(0)=0，z(0)=0.01。t∈[0，80]，微分方程右端函數: