一階常微分方程歐拉法與梯形公式局部截?cái)嗾`差與p階精度Range

一階常微分方程歐拉法與梯形公式局部截?cái)嗾`差與p階精度Range-Kutta公式常微分方程MATLAB求解《數(shù)值分析》23例1.一階常微分方程求解區(qū)域:0≤x≤1.52/16解曲線的斜率轉(zhuǎn)化為方向余弦例2.Logistic模型

初值條件y(0)=0.2一階常微分方程初值問題:數(shù)值方法——取定離散點(diǎn):

x0<x1<x2<···<xN

其中,

y=y(x)是未知函數(shù),右端函數(shù)f(x,y)是已知函數(shù),初值y0

是已知數(shù)據(jù)。求未知函數(shù)y(x)在離散點(diǎn)處的近似值y1,y2,y3,·····,yN3/164/16取定步長(zhǎng):h,記

xn=x0+nh,(n=1,2,···,N)Euler公式:yn+1=yn

+hf(

xn,

yn

)求近似解:y’=f(x,y)梯形公式:

左矩形公式用數(shù)值積分方法離散化常微分方程6/16預(yù)-校方法又稱為修正的Euler法,算法如下

k1=f(xn

,yn),

k2=f(xn+1,yn+hk1),由梯形公式推出的預(yù)-校方法:7/16設(shè)

yn=y(xn),稱

Rn+1=y(xn+1)-

yn+1為局部截?cái)嗾`差.即由泰勒公式Euler公式:yn+1=

yn+hf

(xn,yn)的局部截?cái)嗾`差y(xn+1)–yn+1=y(xn)–yn+O(h2)=O(h2)9/16Euler公式的局部截?cái)嗾`差記為:

O(h2)稱Euler公式具有1階精度。若局部截?cái)嗾`差為:O(hp

+1)

則稱顯式單步法具有

p階精度。例3證明修正的Euler法具有2階精度14/16將預(yù)測(cè)公式代入得

yn+1=yn

+0.5h[f(xn,yn)+f(xn+1,yn+hf(xn,yn))]

f(xn+1,

yn+hf(xn,yn))=f(xn+h,yn+hf(xn,yn))=f(xn,yn)+h[fx’]n+hf(xn,yn)[fy’]n+O(h2)利用臺(tái)勞展開式

y(xn+1)=y(xn)+hy’(xn)+0.5h2y”(xn)+O(h3)0.5h[f(xn,yn)+f(xn+1,

yn+hf(xn,yn))]=hy’(xn)+0.5h2y”(xn)+O(h3)局部截?cái)嗾`差:y(xn+1)–yn+1=y(xn)–yn=O(h3)故修正的Euler法具有2階精度。11/16yn+1=

yn+

hy’(xn)+0.5h2y”(xn)+O(h3)三階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)12/16例4數(shù)值實(shí)驗(yàn):幾種不同求數(shù)值解公式的誤差比較

n102030

40h0.20.10.0667

0.05RK46.862e-0053.747e-0067.071e-007

2.186e-007RK30.0012

1.529e-0044.517e-0051.906e-005RK20.01230.0026

0.0011

5.9612e-004Euler0.10590.05210.03420.025613/16MATLAB求解常微分方程初值問題命令：

(1)用臨時(shí)函數(shù)定義一階微分方程的右端函數(shù)；(2)用MATLAB命令ode23()求數(shù)值解。使用格式：[T,Y]=ode23('F',Tspan,y0)其中，Tspan=[t0，tN]是常微分方程的求解區(qū)域，y0是解的初值實(shí)驗(yàn)例題1蛇形曲線的常微分方程初值問題

MATLAB數(shù)值求解命令F=inline('1./(1+x.^2)-2*y.^2')；ode23(F,[0,6],0)輸出結(jié)果為圖形

[T,y]=ode23(f,[0,6],0)將得到自變量和函數(shù)的離散數(shù)據(jù)

MATLAB解常微分方程初值問題命令數(shù)值求解命令:[x,y]=ode23('f',[a,b],y0)f=inline('y-x.*y.^2');[x,y]=ode23(f,[0,2],1)符號(hào)求解命令:dsolve('eqn1',...)symsxydsolve('Dy=y-x*y^2','y(0)=1','x')ans=1/(x-1+2*exp(-x))解析解:15/161.創(chuàng)建羅倫茨模型右端函數(shù)的M文件