用Matlab解微分方程.ppt

用Matlab解微分方程,一、微分方程的解析解,求微分方程（組）的解析解用函數(shù)dsolve。,dsolve(方程1,方程2,方程n,初始條件,自變量),結(jié)果：u=tan(t+C1),輸入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x),結(jié)果:y=3*exp(-2*x)*sin(5*x),解,解輸入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t)；,結(jié)果為：,x=exp(2*t)*C1+C2*exp(-t)-C2*exp(2*t)+exp(2*t)*C3-3*exp(-t)y=-C1*exp(-2*t)+exp(2*t)*C1+C2*exp(-2*t)+C2*exp(-t)-C2*exp(2*t)+exp(2*t)*C3-C3*exp(-t)z=-C1*exp(-2*t)+exp(2*t)*C1-C2*exp(2*t)+C2*exp(-2*t)+exp(2*t)*C3,二、微分方程的數(shù)值解,（一）常微分方程數(shù)值解的定義,在生產(chǎn)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解。而在實(shí)際上對(duì)初值問(wèn)題，一般是要求得到解在若干個(gè)點(diǎn)上滿(mǎn)足規(guī)定精確度的近似值，或者得到一個(gè)滿(mǎn)足精確度要求的便于計(jì)算的表達(dá)式。,因此，研究常微分方程的數(shù)值解法是十分必要的。,（二）建立數(shù)值解法的一些途徑,1、用差商代替導(dǎo)數(shù),此即歐拉法。,2、使用數(shù)值積分,此即改進(jìn)的歐拉法。,3、使用泰勒公式,以此方法為基礎(chǔ)，有龍格-庫(kù)塔法、線性多步法等方法。,4、數(shù)值公式的精度,當(dāng)一個(gè)數(shù)值公式的截?cái)嗾`差可表示為O（hk+1）時(shí)（k為正整數(shù)，h為步長(zhǎng)），稱(chēng)它是一個(gè)k階公式。,k越大，則數(shù)值公式的精度越高。,歐拉法是一階公式，改進(jìn)的歐拉法是二階公式。龍格-庫(kù)塔法有二階公式和四階公式。線性多步法有四階阿達(dá)姆斯外插公式和內(nèi)插公式。,（三）用Matlab軟件求常微分方程的數(shù)值解,t，x=solver（f,ts,x0,options）,1、在解n個(gè)未知函數(shù)的方程組時(shí)，x0和x均為n維向量，m-文件中的待解方程組應(yīng)以x的分量形式寫(xiě)成.,2、使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組.,注意:,1、建立m-文件vdp1000.m如下：functiondy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)2)*y(2)-y(1);,2、取t0=0，tf=3000，輸入命令：T,Y=ode15s(vdp1000,03000,20);plot(T,Y(:,1),-),3、結(jié)果如圖,解1、建立m-文件rigid.m如下：functiondy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);,2、取t0=0，tf=12，輸入命令：T,Y=ode45(rigid,012,011);plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+),3、結(jié)果如圖,圖中，y1的圖形為實(shí)線，y2的圖形為“*”線，y3的圖形為“+”線.,導(dǎo)彈追蹤問(wèn)題,設(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)彈將它擊中？,解法一（解析法）,軌跡圖見(jiàn)程序chase1,解法二(數(shù)值解),1.建立m-文件eq1.mfunctiondy=eq1(x,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1/5*sqrt(1+y(2)2)/(1-x);,2.取x0=0，xf=0.9999，建立主程序ff6.m如下：x0=0，xf=0.9999x,y=ode15s(eq1,x0 xf,00);plot(x,y(:,1),b.),結(jié)論:導(dǎo)彈大致在（1，0.21）處擊中乙艦,令y1=y,y2=y1，將方程（3）化為一階微分方程組。,解法三(建立參數(shù)方程求數(shù)值解),設(shè)時(shí)刻t乙艦的坐標(biāo)為(X(t)，Y(t)，導(dǎo)彈的坐標(biāo)為(x(t)，y(t).,3因乙艦以速度v0沿直線x=1運(yùn)動(dòng)，設(shè)v0=1，則w=5，X=1，Y=t,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)2+(t-y(2)2);dy(2)=5*(t-y(2)/sqrt(1-y(1)2+(t-y(2)2);,取t0=0，tf=2，建立主程序chase2.m如下：t,y=ode45(eq2,02,00);Y=0:0.01:2;plot(1,Y,-)，holdonplot(y(:,1),y(:,2),*),5.結(jié)果見(jiàn)圖1,導(dǎo)彈大致在（1，0.2）處擊中乙艦，與前面的結(jié)論一致.,圖1,圖2,在chase2.m中，按二分法逐步修改tf，即分別取tf=1，0.5,0.25,直到tf=0.21時(shí)，得圖2.,結(jié)論：時(shí)刻t=0.21時(shí)，導(dǎo)彈在（1，0.21）處擊中乙艦。,慢跑者與狗,一個(gè)慢跑者在平面上沿橢圓以恒定的速率v=1跑步,設(shè)橢圓方程為:x=10+20cost,y=20+5sint.突然有一只狗攻擊他.這只狗從原點(diǎn)出發(fā),以恒定速率w跑向慢跑者,狗的運(yùn)動(dòng)方向始終指向慢跑者.分別求出w=20,w=5時(shí)狗的運(yùn)動(dòng)軌跡.,1.模型建立,設(shè)時(shí)刻t慢跑者的坐標(biāo)為(X(t)，Y(t)，狗的坐標(biāo)為(x(t)，y(t).,則X=10+20cost,Y=20+15sint,狗從(0,0)出發(fā),與導(dǎo)彈追蹤問(wèn)題類(lèi)似，建立狗的運(yùn)動(dòng)軌跡的參數(shù)方程:,2.模型求解,(1)w=20時(shí),建立m-文件eq3.m如下:functiondy=eq3(t,y)dy=zeros(2,1);dy(1)=20*(10+20*cos(t)-y(1)/sqrt(10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);dy(2)=20*(20+15*sin(t)-y(2)/sqrt(10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);,取t0=0，tf=10，建立主程序chase3.m如下：t0=0;tf=10;t,y=ode45(eq3,t0tf,00);T=0:0.1:2*pi;X=10+20*cos(T);Y=20+15*sin(T);plot(X,Y,-)holdonplot(y(:,1),y(:,2),*),在chase3.m，不斷修改tf的值,分別取tf=5,2.5,3.5,至3.15時(shí),狗剛好追上慢跑者.,建立m-文件eq4.m如下:functiondy=eq4(t,y)dy=zeros(2,1);dy(1)=5*(10+20*cos(t)-y(1)/sqrt(10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);dy(2)=5*(20+15*sin(t)-y(2)/sqrt(10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);,取t0=0，tf=10，建立主程序chase4.m如下：t0=0;tf=10;t,y=ode45(eq4,t0tf,00);T=0:0.1:2*pi;X=10+20*cos(T);Y=20+15*sin(T);plot(X,Y,-)holdonplot(y(:,1),y(:,2),*),在chase3.m，不斷修改tf的值,分別取tf=20,40,80,可以看出,狗永遠(yuǎn)追不上慢跑者.,(2)w=5時(shí),地中海鯊魚(yú)問(wèn)題,意大利生物學(xué)家Ancona曾致力于魚(yú)類(lèi)種群相互制約關(guān)系的研究，他從第一次世界大戰(zhàn)期間,地中海各港口捕獲的幾種魚(yú)類(lèi)捕獲量百分比的資料中，發(fā)現(xiàn)鯊魚(yú)等的比例有明顯增加（見(jiàn)下表），而供其捕食的食用魚(yú)的百分比卻明顯下降.顯然戰(zhàn)爭(zhēng)使捕魚(yú)量下降，食用魚(yú)增加，鯊魚(yú)等也隨之增加，但為何鯊魚(yú)的比例大幅增加呢？,他無(wú)法解釋這個(gè)現(xiàn)象，于是求助于著名的意大利數(shù)學(xué)家V.Volterra，希望建立一個(gè)食餌捕食系統(tǒng)的數(shù)學(xué)模型，定量地回答這個(gè)問(wèn)題.,該模型反映了在沒(méi)有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)系，沒(méi)有考慮食餌和捕食者自身的阻滯作用，是Volterra提出的最簡(jiǎn)單的模型.,首先，建立m-文件shier.m如下：functiondx=shier(t,x)dx=zeros(2,1);dx(1)=x(1)*(1-0.1*x(2);dx(2)=x(2)*(-0.5+0.02*x(1);,其次，建立主程序shark.m如下：t,x=ode45(shier,015,252);plot(t,x(:,1),-,t,x(:,2),*)plot(x(:,1),x(:,2),求解結(jié)果：,左圖反映了x1（t）與x2（t）的關(guān)系?？梢圆聹y(cè)：x1（t）與x2（t）都是周期函數(shù)。,模型（二）考慮人工捕獲,設(shè)表示捕獲能力的系數(shù)為e，相當(dāng)于食餌的自然增長(zhǎng)率由r1降為r1-e，捕食者的死亡率由r2增為r2+e,模型求解:,1、分別用m-文件shier1.m和shier2.m定義上述兩個(gè)方程,2、建立主程序shark1.m,求解兩個(gè)方程，并畫(huà)出兩種情況下鯊魚(yú)數(shù)在魚(yú)類(lèi)總數(shù)中所占比例x2(t)/x1(t)+x2(t),實(shí)線為戰(zhàn)前的鯊魚(yú)比例，“*”線為戰(zhàn)爭(zhēng)中的鯊魚(yú)比例,結(jié)論：戰(zhàn)爭(zhēng)中鯊魚(yú)的比例比戰(zhàn)前高！,實(shí)驗(yàn)作業(yè),1.一個(gè)小孩借助長(zhǎng)度為a的硬棒，拉或推某玩具.此小孩沿某曲線行走，計(jì)算并畫(huà)出玩具的軌跡.,2.討論資金積累、國(guó)民收入、與人口增長(zhǎng)的關(guān)系.（1）若國(guó)民平均收入x與按人口平均資金積累y成正比，說(shuō)明僅當(dāng)總資金積累的相對(duì)增長(zhǎng)率k大于人口的相對(duì)增長(zhǎng)率r時(shí)，國(guó)民平均收入才是增長(zhǎng)的.（2）作出k(x)和r(x)的示意圖，分析人口激增會(huì)引起什么后果.,