一维对流方程在三种差分格式下的解

一维对流方程在三种差分格式下的数值解一、实验目的用数值方法计算一维对流方程在A、B、C三种差分格式下的解。∆x取为0.05，∆𝑥∆𝑡取值为0.5,1,2.并作相关讨论𝜕𝜁𝜕𝑡+𝜕𝜁𝜕𝑥=0,−8≤𝑥≤8,𝑡0𝜁(𝑥,0)=1+𝑥,−1≤𝑥≤0𝜁(𝑥,0)=1−𝑥,0≤𝑥≤1𝜁(𝑥,0)=0,𝑥−1或𝑥1二、实验要求1.学会对MS-FORTRAN的基本操作2.用Fortran编写程序计算一维对流方程在A、B、C三种格式下的解。3.讨论各种格式下不同的∆𝑡∆𝑥值的差分格式解的特点。三、实验步骤1.写出一维对流方程的三种格式的差分方程：A格式：𝜁𝑖𝑛+1=𝜁𝑖𝑛+∆𝑡2∆𝑥(𝜁𝑖+1𝑛−𝜁𝑖−1𝑛)B格式：𝜁𝑖𝑛+1=𝜁𝑖𝑛+∆𝑡∆𝑥(𝜁𝑖+1𝑛−𝜁𝑖𝑛)C格式：𝜁𝑖𝑛+1=𝜁𝑖𝑛+∆𝑡∆𝑥(𝜁𝑖𝑛−𝜁𝑖−1𝑛)2.根据上述差分方程，写出求解的fortran程序（附件1,2）;3.根据fortran程序求出三种格式下的解;4.比较各种格式下不同的∆𝑡∆𝑥的数值解，讨论它们的特点。四、实验结果分析利用上述fortran程序的计算结果，可以观察计算数值，初步判断差分格式的收敛性和稳定性；然后利用EXCEL或MATLAB进行绘图，更形象地体现计算结果。以下为利用计算数据在MATLAB所绘制的图（为了使图更清晰，对比更明显，只画了x∈[-3,3]的图）：11.∆𝑡∆𝑥=0.5(∆𝑥∆𝑡=2)A格式：B格式：2C格式：从上述三图看出，A格式收敛不稳定，B格式不收敛，C格式收敛稳定，但有数据耗散的情况。总体来说，在∆𝑡∆𝑥=0.5的情况下，C格式计算精度最高。2.∆𝑡∆𝑥=1(∆𝑥∆𝑡=1)A格式：3B格式：C格式：从上述三图可以看出，A格式和B格式依然不收敛，C格式稳定收敛，没有发生数值耗散，解出的结果最精确。43.∆𝑡∆𝑥=2(∆𝑥∆𝑡=0.5)A格式：B格式：5C格式：从上述三图可以看出，A、B、C格式均不收敛。因为∆𝑡∆𝑥的值过大，导致时间网格过大，所以三种格式都不收敛。五、实验结论在∆𝑡∆𝑥为0.5,1时，A格式收敛不稳定；在∆𝑡∆𝑥为2时，A不收敛。在∆𝑡∆𝑥为0.5,1,2时，B格式均不收敛。在而在∆𝑡∆𝑥=0.5时，C格式收敛，但会有数据耗散；在∆𝑡∆𝑥=1时，C格式稳定收敛；但在∆𝑡∆𝑥=2时，C格式不收敛。该结论与收敛性分析和稳定性分析结果一致。附件：（1）Main.f90文件!使用说明!注意：该程序必须与data.txt在同一文件夹内!若要修改gamma和deltax，可在data.txt中分别修改GAMMA和DX的值!计算完成后请在该程序项目所在文件夹下寻找Output_FormA(B或C).txt文件查看计算结果PROGRAMMAIN!声明变量!TN时间网格数，XN空间网格数!DX，DT/DX(GAMMA)!RB右边界，LB左边界IMPLICITNONEINTEGER::TN,XN,I,JREAL::DX,GAMMAREAL::LB,RBREAL,ALLOCATABLE::U(:,:),X(:)6LOGICALALIVECHARACTER(len=79)::TEMPCHARACTER(len=2)::FORM!声明全局变量commonTN,XN,DX,GAMMA,RB,LBWRITE(*,*)*************Computing***********!打开数据文件data.txt录入6个参数INQUIRE(FILE='data.txt',EXIST=ALIVE)IF(ALIVE)THENOPEN(UNIT=10,FILE=data.txt,STATUS='OLD',ACTION='READ')READ(10,(a))TEMP!标题信息READ(10,*)TN,DX,GAMMA,LB,RB,FORM!读取6个参数WRITE(*,*)ReadsuccessfullyELSEWRITE(*,*)Thedocumentdoesn'texist.STOPENDIF!计算空间网格数XN=(RB-LB)/DX+1!确定U矩阵和X向量的大小ALLOCATE(U(XN,TN+1))ALLOCATE(X(XN))!定义U初值为0U=0!计算t=0时，U的函数值U（:,1）DOI=1,XNX(I)=LB+(I-1)*DXIF(X(I)=-1.AND.X(I)=0)THENU(I,1)=1+X(I)ELSEIF(X(I)0.AND.X(I)=1)THENU(I,1)=1-X(I)ELSEU(I,1)=0ENDIFENDDO!调用子程序计算求解IF(FORM==A.or.FORM==a)THENCALLFormA(U,X)ELSEIF(FORM==B.or.FORM==b)THENCALLFormB(U,X)ELSEIF(FORM==C.or.FORM==c)THENCALLFormC(U,X)ELSEWRITE(*,*)ERROR:Formfault7stopENDIFstopENDPROGRAM!***************子程序A********************SubroutineFormA(U,X)implicitnoneINTEGER::TN,XN,I,JREAL::DX,GAMMAREAL::LB,RBreal::U(XN,TN),X(XN)commonTN,XN,DX,GAMMA,RB,LB!打开输出文件用于输出数据OPEN(UNIT=9,FILE=Output_FormA.txt)!计算T=TN的值WRITE(9,*)差分格式：格式AWRITE(9,*)DT/DX=,GAMMADOJ=2,TNWRITE(9,*)TIME=,J-1WRITE(9,(A5,10X,A1))坐标X,UDOI=1,XNIF(I==1)THEN!边界点置为相邻点的值U(I,J)=U(I+1,J-1)ELSEIF(I==XN)thenU(I,J)=U(I-1,J-1)ELSEU(I,J)=U(I,J-1)-GAMMA*(U(I+1,J-1)-U(I-1,J-1))/2ENDIFWRITE(9,*)X(I),U(I,J)ENDDOENDDOWRITE(*,*)***********Finish**********write(*,*)请在程序文件夹下打开Output_FormA.txt查看结果returnendsubroutine!***************子程序B********************SubroutineFormB(U,X)implicitnoneINTEGER::TN,XN,I,JREAL::DX,GAMMAREAL::LB,RB8real::U(XN,TN),X(XN)commonTN,XN,DX,GAMMA,RB,LB!打开输出文件用于输出数据OPEN(UNIT=9,FILE=Output_FormB.TXT)!计算T=TN的值WRITE(9,*)差分格式：格式BWRITE(9,*)DT/DX=,GAMMADOJ=2,TNWRITE(9,*)TIME=,J-1WRITE(9,(A5,10X,A1))坐标X,UDOI=1,XNIF(I==1)THEN!边界点置为相邻点的值U(I,J)=U(I+1,J-1)ELSEIF(I==XN)thenU(I,J)=U(I-1,J-1)ELSEU(I,J)=U(I,J-1)-GAMMA*(U(I+1,J-1)-U(I,J-1))ENDIFWRITE(9,(F6.3,6X,F6.3))X(I),U(I,J)ENDDOENDDOWRITE(*,*)***********Finish**********write(*,*)请在程序文件夹下打开Output_FormB.txt查看结果returnendsubroutine!***************子程序C********************SubroutineFormC(U,X)implicitnoneINTEGER::TN,XN,I,JREAL::DX,GAMMAREAL::LB,RBreal::U(XN,TN),X(XN)commonTN,XN,DX,GAMMA,RB,LB!打开输出文件用于输出数据OPEN(UNIT=9,FILE=Output_FormC.txt)!计算T=TN的值WRITE(9,*)差分格式：格式CWRITE(9,*)DT/DX=,GAMMADOJ=2,TNWRITE(9,*)TIME=,J-1WRITE(9,(A5,10X,A1))坐标X,UDOI=1,XNIF(I==1)THEN!边界点置为相邻点的值9U(I,J)=U(I+1,J-1)ELSEIF(I==XN)thenU(I,J)=U(I-1,J-1)ELSEU(I,J)=U(I,J-1)-GAMMA*(U(I,J-1)-U(I-1,J-1))ENDIFWRITE(9,(F6.3,6X,F6.3))X(I),U(I,J)ENDDOwrite(9,(80X))ENDDOWRITE(*,*)***********Finish**********write(*,*)请在程序文件夹下打开Output_FormC.txt查看结果returnendsubroutine(2)data.txtTNDXGAMMALBRBFORM50.050.5-88A附件说明：主要更改TN,GAMMA和FORM的值。分别为时间网格数，DT_DX,差分格式其中，TN取任何大于1的整数；GAMMA可取值为0.5,1,2;FORM可取为A,B,C