資源描述:
《拋物型方程數(shù)值解.pdf》由會員上傳分享,免費在線閱讀,更多相關內容在行業(yè)資料-天天文庫。
1、數(shù)學與計算科學學院實驗報告實驗項目名稱拋物型微分方程數(shù)值解所屬課程名稱微分方程數(shù)值解實驗類型驗證性實驗日期2015-4-23班級信計12-2班學號201253100212姓名黃全林成績一、實驗概述:【實驗目的】掌握拋物型方程的有限差分法,并學會應用有限差分法求解拋物型方程數(shù)值解的MATLAB實現(xiàn)。本文利用向前差分法求解拋物型方程的數(shù)值解?!緦嶒炘怼繏佄镄头匠痰南惹安罘址ㄇ蠼猓悍匠屉x散:n?1nnnnu?uu?2u?ujjj?1jj?1??sint,j?1,???,J?1,n?1,???,N2n?h邊值條件:n?1n
2、nnnu0?u0u1?2u0?u?1nn??sin,tu?u,j?0,2n?11?hn?1nnnnuJ?uJuJ?1?2uJ?uJ?1nnj?J,?2?sin,tunJ?1?uJ?1.?h初值條件:u?xj,0??cos?xj,j?0,1,???,J.【實驗環(huán)境】1.硬件環(huán)境第1頁共17頁2.2.軟件環(huán)境MATLAB7.0二、實驗內容:【實驗過程】(實驗步驟)實驗任務求解一維拋物方程的初邊值問題:2?u?u??sin,0t?x?1,t?0,2?x?xux?0,t??ux?1,t??0,t?0,u?x,0??cos?x
3、,0?x?1.2??t?精確解:u?ecos?x??1cos.?t??利用向前差分法求解利用MATLAB進行求解,編輯函數(shù)文件hql_xiangqianchafen.m。源程序見附錄。編輯調用函數(shù)hql_xiangqianchafen(h,m,n,kmax,ep)的腳本文件,并作出相應的求解曲面、精確解曲面和誤差曲面圖形。hql_paowufangcheng.m源程序見附錄。運行上述hql_paowufangcheng.m文件:(1)當空間步長為0.1,時間步長為0.005時,結果如下:第2頁共17頁求解曲面:精確解
4、曲面:誤差曲面:第3頁共17頁通過相應曲面的對比,可以發(fā)現(xiàn),此時精確解與數(shù)值解圖形逼近效果明顯趨同。下面分別固定空間x和時間t,觀察數(shù)值解與精確解的逼近程度。A.時間固定,空間x與u的關系:數(shù)值解:精確解:第4頁共17頁誤差:數(shù)值解與精確解在x方向上的走勢趨同。B.空間x固定,觀察時間t與u的走勢,結果如下:第5頁共17頁數(shù)值解:精確解:誤差:第6頁共17頁可以看出在時間方向基本趨同。(2)縮小步長,觀察對輸出的影響,取空間步長為0.1,時間步長為0.001,輸出結果為:求解曲面:第7頁共17頁精確解曲面:誤差曲面:
5、第8頁共17頁通過相應曲面的對比,可以發(fā)現(xiàn),此時精確解與數(shù)值解圖形逼近效果明顯趨同。下面分別固定空間x和時間t,觀察數(shù)值解與精確解的逼近程度。B.時間固定,空間x與u的關系:數(shù)值解:第9頁共17頁精確解:誤差:數(shù)值解與精確解在x方向上的走勢趨同。第10頁共17頁B.空間x固定,觀察時間t與u的走勢,結果如下:數(shù)值解:精確解:第11頁共17頁誤差:可以看出在時間方向基本趨同。第12頁共17頁【實驗小結】(收獲體會)通過此次上機親自體驗,使我對于數(shù)值計算有了更加深刻的認識。同時,我也感覺自己的編程能力還有待提高,這對于數(shù)
6、值計算是至關重要。這直接影響到了程序的效率以及程序的運行時間以及計算精度。在以后的學習中,更應該加強這方面的鍛煉。三、指導教師評語及成績:評語等級評語及優(yōu)良中不及格格1.實驗報告按時完成,字跡清楚,文字敘述流暢,邏輯性強2.實驗方案設計合理3.實驗過程(實驗步驟詳細,記錄完整,數(shù)據合理,分析透徹)4實驗結論正確.成績:指導教師簽名:批閱日期:附錄:源程序第13頁共17頁函數(shù)文件:文件名:hql_xiangqianchafen.mfunction[puext]=hql_xiangqianchafen(h1,h2,m,n
7、)%解拋物線型一維方程向前歐拉格式(Ut-aUxx=f(x,t),a>0)%不用解線性方程組,由下一層(時間層)的值就直接得到上一層的值%m,n為x,t方向的網格數(shù),例如(2-0)/0.01=200;%e為誤差,p為精確解u=zeros(n+1,m+1);x=0+(0:m)*h1;t=0+(0:n)*h2;for(i=1:n+1)u(i,1)=0;u(i,m+1)=0;endfor(i=1:m+1)u(1,i)=cos(pi*x(i));endfor(i=1:n+1)for(j=1:m+1)f(i,j)=0;ende
8、ndr=h2/(h1*h1);%此處r=a*h2/(h1*h1);a=1要求r<=1/2差分格式才穩(wěn)定for(i=1:n)for(j=2:m)u(i+1,j)=(1-2*r)*u(i,j)+r*(u(i,j-1)+u(i,j+1))+h2*f(i,j);endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(-pi*p