一維對(duì)流方程在A、B、C三種差分格式

1、、上機(jī)目的用數(shù)值方法計(jì)算一維對(duì)流方程在 A B C三種差分格式下的解。X取為0.05.丄取值為0.5,t1, 2。并作相關(guān)討論。、實(shí)驗(yàn)原理t(X, 0)(X, 0)(X, 0)-10, -8 X1 X三、上機(jī)要求：1. 學(xué)會(huì)對(duì)MS-FORTRA的基本操作。2. 用Fortran編寫程序計(jì)算一維對(duì)流方程在 A B C三種格式下的解。3. 討論各種格式下不同的 亠 值的差分格式解的特點(diǎn)。X四、實(shí)驗(yàn)程序以A格式為例，對(duì)微分方程進(jìn)行離散化，得出A格式的差分解的表達(dá)式:B、C格式同理可以寫出。由此編寫如下的Fortran程序。注：除了循環(huán)時(shí)間層的計(jì)算公式略有不同外，三個(gè)程序沒有區(qū)別，因此這里只用一個(gè)主程

2、序，并根據(jù) 格式選擇的不同，采用條件語句判斷執(zhí)行哪一部分循環(huán)體。!空間節(jié)點(diǎn)321, dx=0.05_輸岀依次為(時(shí)間，空間，數(shù)值)_1program ma in Iimp licit none |real dx_dt!定義剛加的值integer abc,r t,i,j,k!定義變量，abc為格式類型，r t為時(shí)間網(wǎng)格數(shù)，其余為循環(huán)變量real , allocatable :s(:,:)!定義存儲(chǔ)矩陣 sIwrite (*,*)“ 輸入dx dt=0.5,1,2Iread (*,*) dx_dtIwrite (*,*)選擇格式，A,B,C分別輸入1,2或3Iread (*,*) abc!根據(jù)格式

3、選擇生成相應(yīng)的文件Iif (abc=1) thenIopen(un it=8,file= out_a.csv)Ielseif (abc=2) thenIopen(un it=8,file=out_b.csv ) |Ielseif (abc=3) then |Iopen(un it=8,file= out_c.csv )I endifr t=160/dx dt !計(jì)算時(shí)間網(wǎng)格總數(shù)allocate (s(r t+1,321)!分配存儲(chǔ)矩陣的空間精選文庫!第一層賦初值do i=1,140,1s(1,i)=0write (8,*) 1 , ,i, , ,s(1,i)end do Ido i=141,1

4、61,1s(1,i)=1+0.05*(i-161) write (8,*) 1 , ,i, , ,s(1,i)end dodo i=162,181,1s(1,i)=1-0.05*(i-161)write (8,*) 1，“，“，i, , ,s(1,i)|end dodo i=182,321,1s(1,i)=0write (8,*) 1 , ,i, , ,s(1,i)end do!循環(huán)時(shí)間層，根據(jù)格式的選擇來判斷執(zhí)行哪一部分if (abc=1) thendo i=2,r_t,1do j=i,322-i,1s(i,j)=s(i-1,j)-(s(i-1,j+1)-s(i-1,j-1)/(dx dt*

5、2)write (8,*)i, ,j,1end do |,s(i,j)do k=1,i-1,1!余下部分賦值0,下同s(i,k)=0write (8,*)i, ,k,s(i,k)end do1do k=322-i,321,1|s(i,k)=0write (8,*)i, ,k,s(i,k)end doend doelseif (abc=2) thendo i=2,r_t+1,1do j=1,322-i,11s(i,j)=s(i-1,j)-(s(i-1,j+1)-s(i-1,j)/dx_dtwrite (8,*)i,j, , ,s(i,j)end dodo k=322-i,321,11s(i,k)

6、=01write (8,*)i,k, , ,s(i,k)end doend do | elseif (abc=3) thendo i=2,r t+1,1do j=i,321,1s(i,j)=s(i-1,j)-(s(i-1,j)-s(i-1,j-1)/dx dtwrite (8,*)i, ,j, , ,s(i,j)end do do k=1,i-1,1s(i,k)=0write (8,*)i, ,k, , ,s(i,k)|end do end do endif |!完成提示write （*,*）數(shù)據(jù)已輸岀至源目錄P ause sto p end pr ogram五、實(shí)驗(yàn)結(jié)果及分析程序運(yùn)行后在對(duì)應(yīng)

7、目錄下生成csv表格文件，根據(jù)輸入的F的值不同生成對(duì)應(yīng)的網(wǎng)格并計(jì)算各節(jié)Ar點(diǎn)數(shù)值。這里以2=1為例，就A、B、C三種格式的結(jié)果進(jìn)行分析。首先我們用MATLAB件畫出初始波形，以便于對(duì)比（這里三個(gè)格式?jīng)]有區(qū)別）A格式，根據(jù)輸出結(jié)果, 每隔10單位取一點(diǎn)（即取選取不同時(shí)間節(jié)點(diǎn)用MATLAB軟件繪制波形圖像，這里從初始時(shí)刻開始t=1,11,21， 畫在同一張圖里，可明顯看出隨時(shí)間變化的趨勢(shì)）。波形畫出如下：（下一頁第一張）從圖中看出，隨著時(shí)間推移A格式呈現(xiàn)明顯的發(fā)散狀態(tài)。B格式圖像如下（下一頁第二張）：從圖中看出，隨著時(shí)間推移B格式也呈現(xiàn)明顯的發(fā)散狀態(tài)。1.50.5-11.5-27C格式圖像如下:0-810.90.00.70.60 50.40.30.2-6-4-202468對(duì)比之下明顯可以看出，隨著時(shí)間的推移， C格式能較好的保證收斂性。同樣的，附上冷=0占以及= 3時(shí)的作圖結(jié)果，從圖中可以看出 三=砧 時(shí)三種格式都呈發(fā)散狀給收斂性，A、B 發(fā)散狀態(tài)。的波形圖B、C格式， 序數(shù)0-321位詈=3時(shí)C格式A精選文庫11-4F面是詈=2