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

解一階常微分方程歐拉法Range-Kutta公式一階微分方程組與二階方程常微分方程邊值問(wèn)題線性多步法簡(jiǎn)介《數(shù)值分析》第八章解一階常微分方程歐拉法《數(shù)值分析》第八章一階常微分方程初值問(wèn)題:數(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ù)值逼近.一階常微分方程初值問(wèn)題:數(shù)值方法——取離散點(diǎn):x0<例1.常微分方程與向量場(chǎng)平面區(qū)域:0≤x≤1.5,0≤y≤1.5,斜率方向余弦例1.常微分方程與向量場(chǎng)平面區(qū)域:0≤x≤1.5,取步長(zhǎng)h,記點(diǎn):在第n個(gè)點(diǎn)處一階向前差商:歐拉公式初值問(wèn)題:取步長(zhǎng)h,記點(diǎn):在第n個(gè)點(diǎn)處一階向前差商:例2.Logistic模型

解析解:歐拉公式例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)取點(diǎn)(xn,yn),(n=0,1,2y’=f(x,y)梯形公式:

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

k2=f(xn+1,yn+hk1),由梯形公式推出的預(yù)-校方法:預(yù)-校方法又稱為修正的Euler法,算法如下由梯形公式推出的設(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階精度。設(shè)yn=y(xn),稱Rn+1=y(xn+1)-若局部截?cái)嗾`差為:O(hp

+1)

則稱顯式單步法具有p階精度。例3證明修正的Euler法具有2階精度證由預(yù)測(cè)公式若局部截?cái)嗾`差為:O(hp+1)則稱顯式單步法具有由Taylor級(jí)數(shù)設(shè)局部截?cái)嗾`差:故修正的Euler法具有2階精度。由Taylor級(jí)數(shù)設(shè)局部截?cái)嗾`差:故修正的Euler法具三階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)三階Range-Kutta公式一般形式y(tǒng)n+1=yn+h例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.0256例4數(shù)值實(shí)驗(yàn):幾種不同求數(shù)值解公式的誤差比較nMATLAB求解常微分方程初值問(wèn)題命令：

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

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ù)

實(shí)驗(yàn)例題1蛇形曲線的常微分方程初值問(wèn)題MATLAB數(shù)值求MATLAB解常微分方程初值問(wèn)題命令數(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))解析解:MATLAB解常微分方程初值問(wèn)題命令數(shù)值求解命令:f=inl兩個(gè)未知函數(shù)的一階常微分方程組常微分方程組的向量形式記兩個(gè)未知函數(shù)的一階常微分方程組常微分方程組的向量形式記歐拉公式:歐拉公式:修改的歐拉公式:(n=0,1,·······,N)經(jīng)典龍格-庫(kù)塔公式:修改的歐拉公式:(n=0,1,·······,N捕食者與被捕食者問(wèn)題

海島上有狐貍和野兔,當(dāng)野兔數(shù)量增多時(shí)，狐貍捕食野兔導(dǎo)致狐群數(shù)量增長(zhǎng);大量兔子被捕食使狐群進(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

捕食者與被捕食者問(wèn)題海島上有狐貍和野兔,當(dāng)野兔數(shù)量增多時(shí)，平面向量場(chǎng):——向量場(chǎng)中過(guò)點(diǎn):(100,20)

的軌線平面向量場(chǎng):——向量場(chǎng)中過(guò)點(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ù)MATLAB命令求解:functionz=fox(t,y)“蝴蝶效應(yīng)”來(lái)源于洛倫茲一次講演。模型如下求微分方程數(shù)值解,

繪出解函數(shù)曲線取=8/3，=10，=28。x(0)=0，y(0)=0，z(0)=0.01。t∈[0，80]，微分方程右端函數(shù):“蝴蝶效應(yīng)”來(lái)源于洛倫茲一次講演。模型如下求微分方程數(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))記向量[y1，y2，y3]=[x，y，z]，創(chuàng)建函數(shù)文分量x

的誤差分量y

的誤差分量z

的誤差分量x的誤差分量y的誤差分量z的誤差二階常微分方程問(wèn)題(簡(jiǎn)諧振動(dòng))(衰減振動(dòng))(受迫振動(dòng))n

階勒讓德方程n階貝塞爾方程二階常微分方程問(wèn)題(簡(jiǎn)諧振動(dòng))(衰減振動(dòng))(受迫振動(dòng))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]);例3.單擺的數(shù)學(xué)模型其中,a=g/L初值條件:[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)模擬程序[t,thata]=ode23('dan',[0,2.755人造衛(wèi)星的軌道模型萬(wàn)有引力定律

地球引力參數(shù):

GM=3.986005×105（km3/s2）

常微分方程人造衛(wèi)星的軌道模型萬(wàn)有引力定律地球引力參數(shù):GM=例4求解邊值問(wèn)題的數(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.高斯消元法例4求解邊值問(wèn)題的數(shù)值方法算例解:取正整數(shù)n,令h=1/–yj-1+(2–h2)yj–yj+1=xjh2

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

y(xn);o

yn–yj-1+(2–h2)yj–yj+1=xj線性多步法(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)/12線性多步法(n=0,1,···)其中,xn+i=x0y’=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)/2y’=f(x,y)在區(qū)間[xn,xn+1]上插值思考題與練習(xí)題1.例舉幾種求解常微分方程初值問(wèn)題的數(shù)值法,簡(jiǎn)述它們各自特點(diǎn)？2.求常微分方程初值問(wèn)題數(shù)值法的局部截?cái)嗾`差與精度有何關(guān)系？3.如何利用一階常微分方程組的數(shù)值方法求解高階常微分方程初值問(wèn)題?4.試列舉出三個(gè)數(shù)學(xué)物理或工程中常見(jiàn)的二階常微分方程(或方程組),簡(jiǎn)述已有的解析方法并與數(shù)值法比較.思考題與練習(xí)題1.例舉幾種求解常微分方程初值問(wèn)題的數(shù)值法,5.將勒讓德第二類橢圓積分轉(zhuǎn)換為常微分方程初值問(wèn)題,并比較常微分方程數(shù)值解法與數(shù)值積分方法相異與相同之處.4.用四階Range-Kutta公式求解下列問(wèn)題有何特殊之處5.將勒讓德第二類橢圓積分轉(zhuǎn)換為常微分方程初值問(wèn)題,并比較解一階常微分方程歐拉法Range-Kutta公式一階微分方程組與二階方程常微分方程邊值問(wèn)題線性多步法簡(jiǎn)介《數(shù)值分析》第八章解一階常微分方程歐拉法《數(shù)值分析》第八章一階常微分方程初值問(wèn)題:數(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ù)值逼近.一階常微分方程初值問(wèn)題:數(shù)值方法——取離散點(diǎn):x0<例1.常微分方程與向量場(chǎng)平面區(qū)域:0≤x≤1.5,0≤y≤1.5,斜率方向余弦例1.常微分方程與向量場(chǎng)平面區(qū)域:0≤x≤1.5,取步長(zhǎng)h,記點(diǎn):在第n個(gè)點(diǎn)處一階向前差商:歐拉公式初值問(wèn)題:取步長(zhǎng)h,記點(diǎn):在第n個(gè)點(diǎn)處一階向前差商:例2.Logistic模型

解析解:歐拉公式例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)取點(diǎn)(xn,yn),(n=0,1,2y’=f(x,y)梯形公式:

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

k2=f(xn+1,yn+hk1),由梯形公式推出的預(yù)-校方法:預(yù)-校方法又稱為修正的Euler法,算法如下由梯形公式推出的設(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階精度。設(shè)yn=y(xn),稱Rn+1=y(xn+1)-若局部截?cái)嗾`差為:O(hp

+1)

則稱顯式單步法具有p階精度。例3證明修正的Euler法具有2階精度證由預(yù)測(cè)公式若局部截?cái)嗾`差為:O(hp+1)則稱顯式單步法具有由Taylor級(jí)數(shù)設(shè)局部截?cái)嗾`差:故修正的Euler法具有2階精度。由Taylor級(jí)數(shù)設(shè)局部截?cái)嗾`差:故修正的Euler法具三階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)三階Range-Kutta公式一般形式y(tǒng)n+1=yn+h例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.0256例4數(shù)值實(shí)驗(yàn):幾種不同求數(shù)值解公式的誤差比較nMATLAB求解常微分方程初值問(wèn)題命令：

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

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ù)

實(shí)驗(yàn)例題1蛇形曲線的常微分方程初值問(wèn)題MATLAB數(shù)值求MATLAB解常微分方程初值問(wèn)題命令數(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))解析解:MATLAB解常微分方程初值問(wèn)題命令數(shù)值求解命令:f=inl兩個(gè)未知函數(shù)的一階常微分方程組常微分方程組的向量形式記兩個(gè)未知函數(shù)的一階常微分方程組常微分方程組的向量形式記歐拉公式:歐拉公式:修改的歐拉公式:(n=0,1,·······,N)經(jīng)典龍格-庫(kù)塔公式:修改的歐拉公式:(n=0,1,·······,N捕食者與被捕食者問(wèn)題

海島上有狐貍和野兔,當(dāng)野兔數(shù)量增多時(shí)，狐貍捕食野兔導(dǎo)致狐群數(shù)量增長(zhǎng);大量兔子被捕食使狐群進(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

捕食者與被捕食者問(wèn)題海島上有狐貍和野兔,當(dāng)野兔數(shù)量增多時(shí)，平面向量場(chǎng):——向量場(chǎng)中過(guò)點(diǎn):(100,20)

的軌線平面向量場(chǎng):——向量場(chǎng)中過(guò)點(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ù)MATLAB命令求解:functionz=fox(t,y)“蝴蝶效應(yīng)”來(lái)源于洛倫茲一次講演。模型如下求微分方程數(shù)值解,

繪出解函數(shù)曲線取=8/3，=10，=28。x(0)=0，y(0)=0，z(0)=0.01。t∈[0，80]，微分方程右端函數(shù):“蝴蝶效應(yīng)”來(lái)源于洛倫茲一次講演。模型如下求微分方程數(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))記向量[y1，y2，y3]=[x，y，z]，創(chuàng)建函數(shù)文分量x

的誤差分量y

的誤差分量z

的誤差分量x的誤差分量y的誤差分量z的誤差二階常微分方程問(wèn)題(簡(jiǎn)諧振動(dòng))(衰減振動(dòng))(受迫振動(dòng))n

階勒讓德方程n階貝塞爾方程二階常微分方程問(wèn)題(簡(jiǎn)諧振動(dòng))(衰減振動(dòng))(受迫振動(dòng))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]);例3.單擺的數(shù)學(xué)模型其中,a=g/L初值條件:[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)模擬程序[t,thata]=ode23('dan',[0,2.755人造衛(wèi)星的軌道模型萬(wàn)有引力定律

地球引力參數(shù):

GM=3.986005×105