一维对流方程在A、B、C三种差分格式

一维对流方程在A、B、C三种差分格式一、上机目的用数值方法计算一维对流方程在A、B、C三种差分格式下的解。x取为0.05.xt取值为0.5，1，2。并作相关讨论。二、实验原理1x1x00)(x,1x010)(x,0x1-x10)(x,0t8x8-,0或xxt三、上机要求：1.学会对MS-FORTRAN的基本操作。2.用Fortran编写程序计算一维对流方程在A、B、C三种格式下的解。3.讨论各种格式下不同的tx值的差分格式解的特点。四、实验程序以A格式为例，对微分方程进行离散化，得出A格式的差分解的表达式：B、C格式同理可以写出。由此编写如下的Fortran程序。注：除了循环时间层的计算公式略有不同外，三个程序没有区别，因此这里只用一个主程序，并根据格式选择的不同，采用条件语句判断执行哪一部分循环体。!空间节点321，dx=0.05输出依次为（时间，空间，数值）programmainimplicitnonerealdx_dt!定义的值integerabc,r_t,i,j,k!定义变量，abc为格式类型，r_t为时间网格数，其余为循环变量real,allocatable::s(:,:)!定义存储矩阵swrite(*,*)输入dx_dt=0.5,1,2read(*,*)dx_dtwrite(*,*)选择格式，A,B,C分别输入1,2或3read(*,*)abc!根据格式选择生成相应的文件if(abc==1)thenopen(unit=8,file='out_a.csv')elseif(abc==2)thenopen(unit=8,file='out_b.csv')elseif(abc==3)thenopen(unit=8,file='out_c.csv')endifr_t=160/dx_dt!计算时间网格总数doj=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)enddodok=1,i-1,1s(i,k)=0write(8,*)i,',',k,',',s(i,k)enddoenddoendif!完成提示write(*,*)'数据已输出至源目录'pausestopendprogram五、实验结果及分析程序运行后在对应目录下生成csv表格文件，根据输入的𝒕的值不同生成对应的网格并计算各节点数值。这里以=𝟏为例，就A、B、C三种格式的结果进行分析。首先我们用MATLAB软件画出初始波形，以便于对比（这里三个格式没有区别）。A格式，根据输出结果，选取不同时间节点用MATLAB软件绘制波形图像，这里从初始时刻开始每隔10单位取一点（即取t=1,11,21，…画在同一张图里，可明显看出随时间变化的趋势）。波形画出如下：（下一页第一张）从图中看出，随着时间推移A格式呈现明显的发散状态。B格式图像如下（下一页第二张）：从图中看出，随着时间推移B格式也呈现明显的发散状态。C格式图像如下：对比之下明显可以看出，随着时间的推移，C格式能较好的保证收敛性。同样的，附上=𝟎.𝟓以及=𝟐时的作图结果，从图中可以看出=𝟎.𝟓时三种格式都呈发散状态，而=𝟐时C格式能保持收敛性，A、B格式都呈发散状态。首先给出=𝟎.𝟓时三种格式的波形图（依次为A、B、C格式，直接以点序数0~321位横坐标）：下面是=𝟐时的波形图：六、实验总结经过分析，可以发现在=𝟎.𝟓时三种格式都不收敛，而=𝟏和=𝟐时C格式收敛，A、B格式发散，这与理论上的结论也相一致。