微分方程數(shù)值解課程設(shè)計(jì)教學(xué)文稿

微分方程數(shù)值解課程

設(shè)計(jì)精品資料精品資料僅供學(xué)習(xí)與交流僅供學(xué)習(xí)與交流,、如有侵權(quán)請(qǐng)聯(lián)系網(wǎng)站刪除謝謝列出第一階段的方程組：設(shè)M0為火箭本身質(zhì)量，m為燃料質(zhì)量，g為重力加速度=9.8m/F,燃料燃燒率為a,空氣阻力的比例系數(shù)為k,F為推進(jìn)力。M0=1400-1080=320kg;V=(F-kvA21F-kv2-4(M。+m)gM0+mg)/(M0+m)M(T+mm=-a ：】一初值土=0引!n=團(tuán)080o由以上各式可以求出t=『時(shí)火箭的速度。再求解第二階段：「 -kv2-MOgV=(-kvA2-M0g)/M0 立= m=0 rh=0可以求出火箭速度降為0的時(shí)刻。將整個(gè)過(guò)程中的時(shí)間向量以及速度向量聯(lián)合起來(lái)，利用第三章所學(xué)插值與數(shù)值積分的方法可以求得任意時(shí)刻火箭的近似高度。2、方法求解常微分方程時(shí)，我分別采用了自己編寫的歐拉公式、改進(jìn)歐拉公式、4級(jí)4階龍格-庫(kù)塔公式，以及MATLA的帶的龍格-庫(kù)塔方法，結(jié)果與分析1、各種公式的對(duì)比首先，我作出了各種不同公式計(jì)算得到的火箭速度隨時(shí)間變化的圖像，圖如下：

從圖中可以看出，各種公式計(jì)算得到的結(jié)果基本一致，為確定其區(qū)別，將圖像放大，放大約2000倍后，得到下圖:從圖中可以看出，自編歐拉公式距離MATLA的帶龍格-庫(kù)塔公式最遠(yuǎn)，精度最差；自編的改進(jìn)歐拉公式和自編的龍格-庫(kù)塔公式結(jié)果基本一致，兩者中自編龍格-庫(kù)塔公式距MATLA的帶龍格-庫(kù)塔公式的結(jié)果稍近。與之前的分析基本一致。然而產(chǎn)生自編龍格-庫(kù)塔公式與MATLA的帶龍格-庫(kù)塔公式之間的差距的原因還未知。由于MATLA的帶龍格-庫(kù)塔公式精度較高，因此以下不在計(jì)算其它幾項(xiàng)公式的結(jié)果。2、第一階段火箭關(guān)閉瞬間的速度v=267261240773261m/s關(guān)閉前瞬間的加速度a=(F-kvA2--kv工-MOgM0g)/M0=0.914286475421320m/sA2 ——__L此時(shí)火箭的高度h=12189.7791507305m3、第二階段由初始速度以及常微分方程可以求得火箭達(dá)到最高點(diǎn)的時(shí)間約為71.3s;此時(shí)火箭的高度.? -加速度a=-9.80000406052111m/sA24、整體過(guò)程下圖為火箭加速度與時(shí)間的關(guān)系。由圖可以看出，火箭一開始的加速度很大，隨著時(shí)間的推移，火箭的燃料有所減少，與此同時(shí)速度有所上升。兩者中前者使火箭的加速度增大，后者使其減小，綜合作用，最終體現(xiàn)為加速度先略有上升，然后慢慢減小。當(dāng)火箭中燃料燃盡時(shí)，火箭喪失推動(dòng)力，因而加速度急劇減小為負(fù)值。此后火箭速度不斷減小，導(dǎo)致火箭所受阻力逐漸減小，因而加速度的絕對(duì)值有所減小，直到最終火箭速度降為零，火箭不受阻力，僅受重力，加速度為重力加速度。

下圖為火箭速度與時(shí)間的關(guān)系:此圖可以看作是由第一幅圖對(duì)時(shí)間積分所得結(jié)果，本圖的斜率對(duì)應(yīng)第一幅圖的值。第一階段火箭加速度為正，因此速度不斷增加，只是增加的速度不斷減慢。燃料燃盡后，加速度變?yōu)樨?fù)值，因此速度開始急劇下降，與此同時(shí)下降的速率不斷減小。最終火箭速度降為00下圖為火箭高度與時(shí)間的關(guān)系:此圖可以看作是由第二幅圖對(duì)時(shí)間積分所得結(jié)果，本圖的斜率對(duì)應(yīng)第一幅圖的值。一開始或減速度較小，因此高度緩慢增加，之后增加的速度不斷提升，直到火箭的燃料耗盡，此后速度不斷減小，但仍為正值，因此火箭繼續(xù)向上飛行，只是高度提升的速度逐漸變慢。知道最終火箭速度降為零。程序清單1、第一階段的微分方程組functiondx=fly(t,y)M0=320;a=18;k=0.4;g=9.8;F=32000;dx=[(-k*y(1)A2+F-(M0+y(2))*g)/(M0+y(2));-a*sign(y(2))];2、第二階段微分方程組functiondx=fly1(t,y)M0=320;k=0.4;g=9.8;F=0;a=0;dx=[(-k*y(1)A2+F-(M0+y(2))*g)/(M0+y(2));-a];3、自編歐拉公式functiony=Euler(fun,ts,x0)NumOfLine=numel(ts);NumOfRow=numel(x0);h=ts(2)-ts(1);y=zeros(NumOfLine,NumOfRow);y(1,：)=x0;fori=1:NumOfLine-1dy=feval(fun,ts(i),y(i,:));y(i+1,:)=y(i,:)+rot90(dy)*h;end4、自編改進(jìn)歐拉公式functiony=ImprovedEuler(fun,ts,x0)NumOfLine=numel(ts);NumOfRow=numel(x0);h=ts(2)-ts(1);y=zeros(NumOfLine,NumOfRow);y(1,：)=x0;fori=1:NumOfLine-1dy=feval(fun,ts(i),y(i,:));y1=y(i,:)+rot90(dy)*h;dy1=feval(fun,ts(i)+h,y1);y(i+1,:)=y(i,:)+h/2*(rot90(dy+dy1));end5、自編龍格-庫(kù)塔公式functiony=myRungeKutta(fun,ts,x0)NumOfLine=numel(ts);NumOfRow=numel(x0);h=ts(2)-ts(1);y=zeros(NumOfLine,NumOfRow);y(1,：)=x0;fori=1:NumOfLine-1k1=feval(fun,ts(i),y(i,:));dy1=y(i,:)+h/2*rot90(k1);k2=feval(fun,ts(i)+h/2,dy1);dy2=y(i,:)+h/2*rot90(k2);k3=feval(fun,ts(i)+h/2,dy2);dy3=y(i,:)+h*rot90(k3);k4=feval(fun,ts(i)+h,dy3);y(i+1,:)=y(i,:)+h/6*rot90(k1+2*k2+2*k3+k4);end6、“自啟動(dòng)”腳本e=0.1;ts1=0:e:60;y00=[0,1080];y01=Euler(@fly,ts1,y00);y11=ImprovedEuler(@fly,ts1,y00);y21=myRungeKutta(@fly,ts1,y00);[t1,y1]=ode45(@fly,ts1,y00);ts2=60:e:80;y00=[y1((60/e+1),1),0];y000=[y01((60/e+1),1),0];y001=[y11((60/e+1),1),0];y002=[y21((60/e+1),1),0];y02=Euler(@fly1,ts2,y000);y12=ImprovedEuler(@fly1,ts1,y001);y22=myRungeKutta(@fly1,ts1,y002);[t2,y2]=ode45(@fly1,ts2,y00);ts=[ts1,ts2];t=[t1;t2];i=1;while(y2(i,1)>0)i=i+1;endy=[y1;y2(1:i-1,:)];y0=[y01;y02(1:i-1,:)];y1=[y11;y12(1:i-1,:)];y2=[y21;y22(1:i-1,:)];j=numel(y)/numel(y00);plot(t(1:j),y(:,1),k,ts(1:j),y0(:,1),'b',ts(1:j),y1(:,1),'y',ts(1:j),y2(:,1),'r');gridon;legend('MATLAB自帶龍格-庫(kù)塔公式','自編的歐拉公式','自編的改進(jìn)歐拉公式','自編的龍格-庫(kù)塔公式');7、自啟動(dòng)腳本2e=0.01;tsi=0:e:60;y0=[0,1080];[t1,y1]=ode45(@fly,ts1,y0);ts2=60:e:80;y0=[y1((60/e+1),1),0];[t2,y2]=ode45(@fly1,ts2,y0);ts=[ts1,ts2];t=[t1;t2];i=1;while(y2(i,1)>0)i=i+1;endy=[yl;y2(l：i-l,：)];y(1,2)=(32000-(320+1080)*9.8)/(320+1080);forj=2:60/e+ih(j)=mySimpson(t(1:j),y(1:j,1));if(j<60/e+1)y(j,2)=(-0.4*y(j,1)A2+32000-(320+y(j,2))*9.8)/(320+y(j,2));elsey(j,2)=(-0.4*y(j,1)A2-320*9.8)/320;endendh=h(1：j);t=t(1：j);plot(t,h,'k');gridon;總結(jié)本文對(duì)常微分方程初值問(wèn)題現(xiàn)有的數(shù)值解法進(jìn)行了綜述研究。主要討論了幾種常用的數(shù)值解法：即歐拉法，改進(jìn)歐拉法，龍格庫(kù)塔方法，阿達(dá)姆斯PECEFMECMK；等。文章最后結(jié)合常見數(shù)值解法，對(duì)較為典型的微分方程模型進(jìn)行數(shù)值求解，探討了上述數(shù)值算法在實(shí)際建模問(wèn)題中的應(yīng)用本次課程設(shè)計(jì)雖然遇到了很多困難和不懂的地方，但是也學(xué)到了很多東西，覺(jué)得還是很有收獲。這次的課程設(shè)計(jì)按照老師的要求借閱了相應(yīng)的教材，然后