用Matlab解微分方程.ppt
《用Matlab解微分方程.ppt》由會員分享,可在線閱讀,更多相關(guān)《用Matlab解微分方程.ppt(28頁珍藏版)》請在裝配圖網(wǎng)上搜索。
用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)和科研中所處理的微分方程往往很復雜且大多得不出一般解。而在實際上對初值問題,一般是要求得到解在若干個點上滿足規(guī)定精確度的近似值,或者得到一個滿足精確度要求的便于計算的表達式。,因此,研究常微分方程的數(shù)值解法是十分必要的。,(二)建立數(shù)值解法的一些途徑,1、用差商代替導數(shù),此即歐拉法。,2、使用數(shù)值積分,此即改進的歐拉法。,3、使用泰勒公式,以此方法為基礎,有龍格-庫塔法、線性多步法等方法。,4、數(shù)值公式的精度,當一個數(shù)值公式的截斷誤差可表示為O(hk+1)時(k為正整數(shù),h為步長),稱它是一個k階公式。,k越大,則數(shù)值公式的精度越高。,歐拉法是一階公式,改進的歐拉法是二階公式。龍格-庫塔法有二階公式和四階公式。線性多步法有四階阿達姆斯外插公式和內(nèi)插公式。,(三)用Matlab軟件求常微分方程的數(shù)值解,t,x=solver(f,ts,x0,options),1、在解n個未知函數(shù)的方程組時,x0和x均為n維向量,m-文件中的待解方程組應以x的分量形式寫成.,2、使用Matlab軟件求數(shù)值解時,高階微分方程必須等價地變換成一階微分方程組.,注意:,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的圖形為實線,y2的圖形為“*”線,y3的圖形為“+”線.,導彈追蹤問題,設位于坐標原點的甲艦向位于x軸上點A(1,0)處的乙艦發(fā)射導彈,導彈頭始終對準乙艦.如果乙艦以最大的速度v0(是常數(shù))沿平行于y軸的直線行駛,導彈的速度是5v0,求導彈運行的曲線方程.又乙艦行駛多遠時,導彈將它擊中?,解法一(解析法),軌跡圖見程序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é)論:導彈大致在(1,0.21)處擊中乙艦,令y1=y,y2=y1,將方程(3)化為一階微分方程組。,解法三(建立參數(shù)方程求數(shù)值解),設時刻t乙艦的坐標為(X(t),Y(t),導彈的坐標為(x(t),y(t).,3因乙艦以速度v0沿直線x=1運動,設v0=1,則w=5,X=1,Y=t,4.解導彈運動軌跡的參數(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é)果見圖1,導彈大致在(1,0.2)處擊中乙艦,與前面的結(jié)論一致.,圖1,圖2,在chase2.m中,按二分法逐步修改tf,即分別取tf=1,0.5,0.25,直到tf=0.21時,得圖2.,結(jié)論:時刻t=0.21時,導彈在(1,0.21)處擊中乙艦。,慢跑者與狗,一個慢跑者在平面上沿橢圓以恒定的速率v=1跑步,設橢圓方程為:x=10+20cost,y=20+5sint.突然有一只狗攻擊他.這只狗從原點出發(fā),以恒定速率w跑向慢跑者,狗的運動方向始終指向慢跑者.分別求出w=20,w=5時狗的運動軌跡.,1.模型建立,設時刻t慢跑者的坐標為(X(t),Y(t),狗的坐標為(x(t),y(t).,則X=10+20cost,Y=20+15sint,狗從(0,0)出發(fā),與導彈追蹤問題類似,建立狗的運動軌跡的參數(shù)方程:,2.模型求解,(1)w=20時,建立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時,狗剛好追上慢跑者.,建立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,可以看出,狗永遠追不上慢跑者.,(2)w=5時,地中海鯊魚問題,意大利生物學家Ancona曾致力于魚類種群相互制約關(guān)系的研究,他從第一次世界大戰(zhàn)期間,地中海各港口捕獲的幾種魚類捕獲量百分比的資料中,發(fā)現(xiàn)鯊魚等的比例有明顯增加(見下表),而供其捕食的食用魚的百分比卻明顯下降.顯然戰(zhàn)爭使捕魚量下降,食用魚增加,鯊魚等也隨之增加,但為何鯊魚的比例大幅增加呢?,他無法解釋這個現(xiàn)象,于是求助于著名的意大利數(shù)學家V.Volterra,希望建立一個食餌捕食系統(tǒng)的數(shù)學模型,定量地回答這個問題.,該模型反映了在沒有人工捕獲的自然環(huán)境中食餌與捕食者之間的制約關(guān)系,沒有考慮食餌和捕食者自身的阻滯作用,是Volterra提出的最簡單的模型.,首先,建立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:x1(t)與x2(t)都是周期函數(shù)。,模型(二)考慮人工捕獲,設表示捕獲能力的系數(shù)為e,相當于食餌的自然增長率由r1降為r1-e,捕食者的死亡率由r2增為r2+e,模型求解:,1、分別用m-文件shier1.m和shier2.m定義上述兩個方程,2、建立主程序shark1.m,求解兩個方程,并畫出兩種情況下鯊魚數(shù)在魚類總數(shù)中所占比例x2(t)/x1(t)+x2(t),實線為戰(zhàn)前的鯊魚比例,“*”線為戰(zhàn)爭中的鯊魚比例,結(jié)論:戰(zhàn)爭中鯊魚的比例比戰(zhàn)前高!,實驗作業(yè),1.一個小孩借助長度為a的硬棒,拉或推某玩具.此小孩沿某曲線行走,計算并畫出玩具的軌跡.,2.討論資金積累、國民收入、與人口增長的關(guān)系.(1)若國民平均收入x與按人口平均資金積累y成正比,說明僅當總資金積累的相對增長率k大于人口的相對增長率r時,國民平均收入才是增長的.(2)作出k(x)和r(x)的示意圖,分析人口激增會引起什么后果.,- 配套講稿:
如PPT文件的首頁顯示word圖標,表示該PPT已包含配套word講稿。雙擊word圖標可打開word文檔。
- 特殊限制:
部分文檔作品中含有的國旗、國徽等圖片,僅作為作品整體效果示例展示,禁止商用。設計者僅對作品中獨創(chuàng)性部分享有著作權(quán)。
- 關(guān) 鍵 詞:
- Matlab 微分方程
裝配圖網(wǎng)所有資源均是用戶自行上傳分享,僅供網(wǎng)友學習交流,未經(jīng)上傳用戶書面授權(quán),請勿作他用。
鏈接地址:http://appdesigncorp.com/p-3447616.html