《數(shù)學(xué)實(shí)驗(yàn) 第4版》課件 4.4 常微分方程的數(shù)值解

第四章數(shù)值分析

實(shí)驗(yàn)4.1插值

實(shí)驗(yàn)4.2離散數(shù)據(jù)的曲線擬合數(shù)學(xué)實(shí)驗(yàn)實(shí)驗(yàn)4.3MATLAB數(shù)值積分與微分實(shí)驗(yàn)4.4常微分方程的數(shù)值解2實(shí)驗(yàn)4.4常微分方程的數(shù)值解一、幾種求常微分方程數(shù)值解的方法二、應(yīng)用舉例實(shí)驗(yàn)4.4常微分方程的數(shù)值解3實(shí)驗(yàn)4.4常微分方程的數(shù)值解實(shí)驗(yàn)?zāi)康耐ㄟ^本實(shí)驗(yàn)了解常微分方程的數(shù)值解的概念，掌握利用MATLAB求常微分方程的數(shù)值解的方法.一、幾種求常微分方程數(shù)值解的方法常微分方程是研究函數(shù)變化規(guī)律的有力工具，但是絕大多數(shù)變系數(shù)方程、非線性方程都是所謂“解不出來”的，于是常微分方程的數(shù)值解法就成為解常微分方程的主要手段．本實(shí)驗(yàn)只考慮初值問題.4實(shí)驗(yàn)4.4常微分方程的數(shù)值解常微分方程初值問題：保證方程的解存在且惟一在一系列離散點(diǎn)上，求的近似值，通常取等步長h，即因此數(shù)值解法得到的近似解是一個(gè)離散的函數(shù)表.5實(shí)驗(yàn)4.4常微分方程的數(shù)值解一、幾種求常微分方程數(shù)值解的方法1.歐拉（Euler）方法歐拉方法是一種最古老、最簡(jiǎn)單而直觀的解微分方程的數(shù)值方法，其基本想法是在小區(qū)間上用差商代替方程左端的導(dǎo)數(shù),而方程右端已知函數(shù)中的x在小區(qū)間上的哪一點(diǎn)取值，則有以下不同的方法：6實(shí)驗(yàn)4.4常微分方程的數(shù)值解（1）向前歐拉公式中的x取小區(qū)間的左端點(diǎn)xn，以分別代替，就得到：向前歐拉公式（2）向后歐拉公式中的x取小區(qū)間的右端點(diǎn)xn+1，就得到：向后歐拉公式上式右端的yn+1未知，故稱為隱式公式，無法用它直接計(jì)算yn+1．7實(shí)驗(yàn)4.4常微分方程的數(shù)值解（3）梯形公式將向前歐拉公式和向后歐拉公式加以平均，得到梯形公式8實(shí)驗(yàn)4.4常微分方程的數(shù)值解（4）改進(jìn)的歐拉公式先由向前歐拉公式算出yn+1的預(yù)測(cè)值，再把它代入梯形公式右端，作為校正，即改進(jìn)的歐拉公式9實(shí)驗(yàn)4.4常微分方程的數(shù)值解它還可寫作上面4個(gè)公式中，我們常用的是便于計(jì)算的向前歐拉公式和改進(jìn)的歐拉公式.10實(shí)驗(yàn)4.4常微分方程的數(shù)值解2.龍格—庫塔（Runge-Kutta）公式在區(qū)間內(nèi)多取幾個(gè)點(diǎn)，就可以構(gòu)造出精度更高的計(jì)算公式，這就是龍格—庫塔公式，它是一類方法的總稱.它是歐拉方法的一種推廣，也是應(yīng)用最廣的求解常微分方程數(shù)值問題的方法.一般的龍格—庫塔方法的形式為p階龍格—庫塔方法其中ai,bij,ci為待定參數(shù)11實(shí)驗(yàn)4.4常微分方程的數(shù)值解（1）二階龍格—庫塔公式其中為待定系數(shù).滿足：上式有4個(gè)未知數(shù)而只有3個(gè)方程，所以解不唯一，存在無窮多個(gè)解，可見二階龍格—庫塔方法是一族公式.

改進(jìn)的歐拉公式12實(shí)驗(yàn)4.4常微分方程的數(shù)值解（2）四階經(jīng)典（標(biāo)準(zhǔn)）龍格—庫塔公式實(shí)際應(yīng)用中，用的最多的是四階龍格—庫塔公式四階龍格—庫塔公式也不止一個(gè)13實(shí)驗(yàn)4.4常微分方程的數(shù)值解3.龍格—庫塔方法的MATLAB實(shí)現(xiàn)在MATLAB中，求微分方程數(shù)值解的函數(shù)調(diào)用格式為[t,x]=solver（’f’,ts,x0,options）自變量值函數(shù)值ode45、ode23、ode113、ode15s、ode23s中的一個(gè)由待解方程寫成的m-文件名=[t0,tf],t0、tf為自變量的初值和終值函數(shù)的初值用于設(shè)定誤差限options=odeset（'reltol',rt,'abstol',at）默認(rèn)相對(duì)誤差為10-3，絕對(duì)誤差為10-6.相對(duì)誤差絕對(duì)誤差14實(shí)驗(yàn)4.4常微分方程的數(shù)值解在諸多的solver函數(shù)中，ode23、ode45都是基于龍格—庫塔公式的函數(shù):ode23是采用的2、3階龍格-庫塔組合算法，ode45是采用的4、5階龍格-庫塔組合算法.ode45比ode23精度更高，尤其成為大部分場(chǎng)合的首選算法.注(1)在解n個(gè)未知函數(shù)的方程組時(shí)，x0和x均為n維向量，m-文件中的待解方程組應(yīng)以x的分量形式寫成.(2)使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組.15例21

求微分方程的數(shù)值解，解先編寫函數(shù)文件：functiondx=fun23(t,x)dx=t+x;再求方程的解實(shí)驗(yàn)4.4常微分方程的數(shù)值解已知精確解為>>ts=0:0.1:1;x0=1;[t,x]=ode45('fun23',ts,x0);y=2*exp(t)-t-1;[t,x,y]plot(t,x,'r-',t,y,'b-.'),gridon;title('SolutionofExample3');xlabel('timet');ylabel('solutionx');legend('x','y');↙ans=01.00001.000016實(shí)驗(yàn)4.4常微分方程的數(shù)值解0.10001.11031.11030.20001.24281.24280.30001.39971.39970.40001.58361.58360.50001.79741.79740.60002.04422.04420.70002.32752.32750.80002.65112.65110.90003.01923.01921.00003.43663.4366

t的值x的值y的值例22

求微分方程組的數(shù)值解.實(shí)驗(yàn)4.4常微分方程的數(shù)值解17解先編寫函數(shù)文件：functiondx=fun24(t,x)dx=[3*t*x(1)+6*x(2);t^2+x(1)*x(2)];再求方程組的解>>ts=0:0.1:1.5;x0=[1,0];[t,x]=ode45('fun24',ts,x0);a=[t,x];plot(t,x(:,1),'r-',t,x(:,2),'b-.'),gridon↙18a=01.000000.10001.01520.00030.20001.06270.00280.30001.14910.00980.40001.28650.02410.50001.49520.04940.60001.80790.09110.70002.27750.15850.80002.99200.27010.90004.10490.4707實(shí)驗(yàn)4.4常微分方程的數(shù)值解1.00005.90410.88741.10009.01972.00691.200015.32486.75821.300037.960170.5533t的值x的值y的值實(shí)驗(yàn)4.4常微分方程的數(shù)值解例23求解描述振蕩器的經(jīng)典的VerderPol微分方程:解原方程等價(jià)為，，.，1920首先編寫函數(shù)文件verderpol.m：functionxprime=verderpol(t,x)globalmu;xprime=[x(2);mu*(1-x(1)^2)*x(2)-x(1)]；再在命令窗口中執(zhí)行：globalmu;mu=7;y0=[1;0];[t,x]=ode45('verderpol',[0,40],y0);x1=x(:,1);x2=x(:,2);plot(t,x1,t,x2)↙實(shí)驗(yàn)4.4常微分方程的數(shù)值解21通過這些繁雜的公式，我們看到數(shù)學(xué)家們?yōu)榱四Ｐ偷耐晟坪透倪M(jìn)付出無數(shù)的心血和努力，這種對(duì)科學(xué)精益求精的探索和創(chuàng)新精神值得我們學(xué)習(xí).同時(shí)，為了更好地應(yīng)用這些成果，我們應(yīng)該具體問題具體分析，選取適當(dāng)?shù)姆椒ń鉀Q問題.實(shí)驗(yàn)4.4常微分方程的數(shù)值解22實(shí)驗(yàn)4.4常微分方程的數(shù)值解二、應(yīng)用舉例導(dǎo)彈追蹤問題設(shè)位于坐標(biāo)原點(diǎn)的甲艦向位于x軸上點(diǎn)A(1,0)處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對(duì)準(zhǔn)乙艦．如果乙艦以最大的速度v0(是常數(shù))沿平行于y軸的直線行駛，導(dǎo)彈的速度是5v0，求導(dǎo)彈運(yùn)行的曲線方程．又乙艦行駛多遠(yuǎn)時(shí)，導(dǎo)彈將它擊中？解法一（解析法）：假設(shè)導(dǎo)彈在t時(shí)刻的位置為乙艦位于.由于導(dǎo)彈頭始終對(duì)準(zhǔn)乙艦，故此時(shí)直線PQ就是導(dǎo)彈的軌跡曲線弧在點(diǎn)P處的切線，如圖所示：23實(shí)驗(yàn)4.4常微分方程的數(shù)值解圖4.14即有即(1)又根據(jù)題意,弧的長度為的5倍，即(2)由(1)，(2)兩式消去t整理得模型:初值條件為:初值問題的解即為導(dǎo)彈的運(yùn)行軌跡:24當(dāng)時(shí),即當(dāng)乙艦航行到點(diǎn)處時(shí)被導(dǎo)彈擊中.實(shí)驗(yàn)4.4常微分方程的數(shù)值解若,則在處被擊中.被擊中時(shí)間為:編寫以下程序可得軌跡圖>>x=0:0.01:1;y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24;plot(x,y,'*')↙25實(shí)驗(yàn)4.4常微分方程的數(shù)值解解法二(數(shù)值解)(1)建立m-文件eq1.mfunctiondy=eq1(x,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);(2)取x0=0，xf=0.9999，求解如下：>>x0=0;xf=0.9999;[x,y]=ode15s('eq1',[x0xf],[00]);plot(x,y(:,1),'b.')holdony=0:0.01:2;plot(1,y,'b*')↙26結(jié)論:導(dǎo)彈大致在（1,0.2）處擊中乙艦.實(shí)驗(yàn)4.4常微分方程的數(shù)值解27實(shí)驗(yàn)4.4常微分方程的數(shù)值解解法三(建立參數(shù)方程求數(shù)值解)設(shè)時(shí)刻乙艦的坐標(biāo)為，導(dǎo)彈的坐標(biāo)為(1)設(shè)導(dǎo)彈速度恒為(2)由于彈頭始終對(duì)準(zhǔn)乙艦，故導(dǎo)彈的速度平行于乙艦與導(dǎo)彈頭位置的差向量,即：28實(shí)驗(yàn)4.4常微分方程的數(shù)值解消去λ得：(3)因乙艦以速度沿直線運(yùn)動(dòng)因此導(dǎo)彈運(yùn)動(dòng)軌跡的參數(shù)方程為：，則設(shè)29實(shí)驗(yàn)4.4常微分方程的數(shù)值解(4)解導(dǎo)彈運(yùn)動(dòng)軌跡的參數(shù)方程建立m-文件eq2.m如下：functiondy=eq2(t,y)dy=zeros(2,1);dy(1)=5*(1-y(1))/sqrt((1-y(1)