一维对流方程在ABC三种差分格式基础教育

1、一、上机目的用数值方法计算一维对流方程在a、b、c三种差分格式下的解。取为0.05. 取值为0.5，1，2。并作相关讨论。二、实验原理三、上机要求：1.学会对ms-fortran的基本操作。 2.用fortran编写程序计算一维对流方程在a、b、c三种格式下的解。3.讨论各种格式下不同的值的差分格式解的特点。四、实验程序以a格式为例，对微分方程进行离散化， 得出 a 格式的差分解的表达式：b、c格式同理可以写出。由此编写如下的fortran程序。注：除了循环时间层的计算公式略有不同外，三个程序没有区别，因此这里只用一个主程序，并根据格式选择的不同，采用条件语句判断执行哪一部分循环体。!空间节点

2、321，dx=0.05 输出依次为（时间，空间，数值）program mainimplicit nonereal dx_dt !定义x/t的值integer abc,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或3 read(*,*) abc !根据格式选择生成相应的文件 if(abc=1) then open(unit=8,file=out_a.csv)

3、 elseif(abc=2) then open(unit=8,file=out_b.csv) elseif(abc=3) then open(unit=8,file=out_c.csv) endifr_t=160/dx_dt !计算时间网格总数allocate(s(r_t+1,321) !分配存储矩阵的空间!第一层赋初值do i=1,140,1s(1,i)=0write(8,*)1,i,s(1,i)end dodo i=141,161,1s(1,i)=1+0.05*(i-161)write(8,*)1,i,s(1,i)end do do i=162,181,1s(1,i)=1-0.05*(i

4、-161)write(8,*)1,i,s(1,i)end do do i=182,321,1s(1,i)=0write(8,*)1,i,s(1,i)end do !循环时间层,根据格式的选择来判断执行哪一部分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*2) write(8,*)i,j,s(i,j) end do do k=1,i-1,1 !余下部分赋值0,下同 s(i,k)=0 write(8,*)i,k,s(i,k) end do do k=322-i,321,1

5、s(i,k)=0 write(8,*)i,k,s(i,k) end do end do elseif(abc=2) then do i=2,r_t+1,1 do j=1,322-i,1 s(i,j)=s(i-1,j)-(s(i-1,j+1)-s(i-1,j)/dx_dt write(8,*)i,j,s(i,j) end do do k=322-i,321,1 s(i,k)=0 write(8,*)i,k,s(i,k) end do end do elseif(abc=3) then do i=2,r_t+1,1 do j=i,321,1 s(i,j)=s(i-1,j)-(s(i-1,j)-s(

6、i-1,j-1)/dx_dt write(8,*)i,j,s(i,j) end do do k=1,i-1,1 s(i,k)=0 write(8,*)i,k,s(i,k) end do end doendif!完成提示write(*,*)数据已输出至源目录pausestopend program五、实验结果及分析程序运行后在对应目录下生成csv表格文件，根据输入的xt的值不同生成对应的网格并计算各节点数值。这里以xt=1为例，就a、b、c三种格式的结果进行分析。首先我们用matlab软件画出初始波形，以便于对比（这里三个格式没有区别）。a格式，根据输出结果，选取不同时间节点用matlab软件绘制波形图像，这里从初始时刻开始每隔10单位取一点（即取t=1,11,21， 画在同一张图里，可明显看出随时间变化的趋势）。波形画出如下：（下一页第一张）从图中看出，随着时间推移a格式呈现明显的发散状态。b格式图像如下（下一页第二张）：从图中看出，随着时间推移b格式也呈现明显的发散状态。c格式图像如下：对比之下明显可以看出，随着时间的推移，c格式能较好的保证收敛性。同样的，附上xt=0.5以及xt=2时的作图结果，从图中可以看出xt=0.5时三种格式都呈发散状态，而xt=2时c格式能保持收敛性，a、b格式都呈发散状态。首先给出xt=0.5时三种格式的波形图（依次为a、b