數(shù)值分析第八章New

解一階常微分方程歐拉法Range-Kutta公式一階微分方程組與二階方程常微分方程邊值問題線性多步法簡介《數(shù)值分析》第八章一階常微分方程初值問題:數(shù)值方法——取離散點(diǎn):

x0<x1<x2<···<xn

······

其中,

y=y(x)是未知函數(shù),f(x,y)是已知函數(shù).初值y0是已知數(shù)據(jù)。求未知函數(shù)y(x1),y(x2),····,y(xn),······的近似值y1,y2,y3,·····,yn·······稱為常微分方程的數(shù)值解。這里是的數(shù)值逼近.例1.常微分方程與向量場平面區(qū)域:0≤x≤1.5,0≤y≤1.5,斜率方向余弦取步長h,記點(diǎn):在第n個(gè)點(diǎn)處一階向前差商:歐拉公式初值問題:例2.Logistic模型

解析解:歐拉公式取點(diǎn)

(xn

,yn),(n=0,1,2,···

)歐拉公式解的幾何解釋:取

x=xn+1得:yn+1=yn+hf(xn,yn)點(diǎn)斜式直線方程:y=yn+(x–xn)f(xn,yn),取點(diǎn)

(xn

,yn)(xn+1,yn+1)y’=f(x,y)梯形公式:

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

k2=f(xn+1,yn+hk1),由梯形公式推出的預(yù)-校方法:設(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)Euler公式的局部截?cái)嗾`差記為:

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

+1)

則稱顯式單步法具有p階精度。例3證明修正的Euler法具有2階精度證由預(yù)測公式由Taylor級數(shù)設(shè)局部截?cái)嗾`差:故修正的Euler法具有2階精度。三階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í)驗(yàn):幾種不同求數(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.0256MATLAB求解常微分方程初值問題命令：

(1)定義一階微分方程的右端函數(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))解析解:兩個(gè)未知函數(shù)的一階常微分方程組常微分方程組的向量形式記歐拉公式:修改的歐拉公式:(n=0,1,·······,N)經(jīng)典龍格-庫塔公式:捕食者與被捕食者問題

海島上有狐貍和野兔,當(dāng)野兔數(shù)量增多時(shí)，狐貍捕食野兔導(dǎo)致狐群數(shù)量增長;大量兔子被捕食使狐群進(jìn)入饑餓狀態(tài)其數(shù)量下降;狐群數(shù)量下降導(dǎo)致兔子被捕食機(jī)會(huì)減少,兔群數(shù)量回升。微分方程模型如下計(jì)算x(t)，y(t)當(dāng)t∈[0，20]時(shí)的數(shù)據(jù)。繪圖并分析捕食者和被捕食者的數(shù)量變化規(guī)律。x(0)=100y(0)=20

平面向量場:——向量場中過點(diǎn):(100,20)

的軌線MATLAB命令求解:Y0=[100,20];[t,Y]=ode23('fox',[0,20],Y0);x=Y(:,1);y=Y(:,2);figure(1),plot(t,x,'b',t,y,'r')figure(2),plot(x,y)functionz=fox(t,y)z(1,:)=y(1)-0.01*y(1).*y(2);z(2,:)=-y(2)+0.02*y(1).*y(2);----y1----y2y1—y2

相位圖定義方程右端函數(shù)“蝴蝶效應(yīng)”來源于洛倫茲一次講演。模型如下求微分方程數(shù)值解,

繪出解函數(shù)曲線取=8/3，=10，=28。x(0)=0，y(0)=0，z(0)=0.01。t∈[0，80]，微分方程右端函數(shù):記向量[y1，y2，y3]=[x，y，z]，創(chuàng)建函數(shù)文件functionz=flo(t,y)z(1,:)=-8*y(1)/3+y(2).*y(3);z(2,:)=-10*(y(2)-y(3));z(3,:)=-y(1).*y(2)+28*y(2)-y(3);用MATLAB命令求解并繪出Y-X平面的投影圖

P0=[0;0;0.01];[T,P]=ode23('flo',[0,80],P0);figure(1),plot(P(:,2),P(:,1))figure(2),comet3(P(:,1),P(:,2),P(:,3))分量x

的誤差分量y

的誤差分量z

的誤差二階常微分方程問題(簡諧振動(dòng))(衰減振動(dòng))(受迫振動(dòng))n

階勒讓德方程n階貝塞爾方程令一階常微分方程組:初值條件:常微分方程組例3.單擺的數(shù)學(xué)模型其中,a=g/L初值條件:(0)=0.4，

’(0)=0第一步:轉(zhuǎn)化為一階方程組令:y1=,y2=

’

初值條件:y1(0)=0.4,y2(0)=0第二步:求解方程組functionf=dan(x,y)f(1,:)=y(2);f(2,:)=-9.8*sin(y(1))/3.2;L=3.2ode23('dan',[0,2],[0.4,0]);[t,thata]=ode23('dan',[0,2.755],[0.6,0]);R=3.2;n=length(t);alpha=thata(:,1);x=R*sin(alpha);y=R*cos(alpha);X=[0,0];Y=[0,-3.5];fork=1:nxk=x(1:k);yk=y(1:k);Xk=x(k);Yk=y(k);plot(xk,-yk,'.-r',Xk,-Yk,'o',[0,Xk],[0,-Yk]),axis([-2.5,2.5,-3.5,0])pause(.5)end單擺的動(dòng)態(tài)模擬程序人造衛(wèi)星的軌道模型萬有引力定律

地球引力參數(shù):

GM=3.986005×105（km3/s2）

常微分方程例4求解邊值問題的數(shù)值方法算例解:取正整數(shù)n,令h=1/(n+1),xj=jh,(j=0,1,···,n+1

).

將常微分方程離散化

整理,得:

–yj-1+(2–h2)yj–yj+1=xjh2

(j=1,2,···,n)y0=0,yn=01.打靶法;2.高斯消元法–yj-1+(2–h2)yj–yj+1=xjh2

(j=1,2,···,n)三對角方程組AY=F—

y(xn);o

yn線性多步法(n=0,1,···)其中,xn+i=x0+(n+i)h,fn+i=f(xn+i,yn+i)局部載斷誤差A(yù)damas顯格式:yn+2=yn+1+h(3fn+1-fn)/2yn+3=yn+2+h(23fn+2-16fn+1+5fn)/12y’=f(x,y)在區(qū)間[xn,xn+1]上插值f(x)≈[(xn+1－x)fn+(x－xn)fn+1]/h二階Adamas顯格式:yn+2=yn+1+h(3fn+1-fn)/2思考題與練習(xí)題1.例舉幾種求