一维对流方程在三种差分格式下的数值解
一、实验目的
用数值方法计算一维对流方程在A 、B 、C 三种差分格式下的解。∆x 取为0.05，∆t为0.5,1,2. 并作相关讨论
ðζðζ+=0, −8≤x≤8 ,t>0 ðtðx
ζ(x,0) =1+x,−1≤x≤0 ζ(x,0) =1−x,0≤x≤1 ζ(x,0) =0,
x1
∆x
二、实验要求
1. 学会对MS-FORTRAN 的基本操作
2. 用Fortran 编写程序计算一维对流方程在A 、B 、C 三种格式下的解。 3. 讨论各种格式下不同的
∆x∆t
三、实验步骤
1. 写出一维对流方程的三种格式的差分方程：
nn
A 格式：ζin+1=ζin+2∆x(ζi+1−ζi−1)
∆t
B 格式：ζin+1=ζin+C 格式：ζin+1=ζin+
∆t∆x∆t∆x
n
(ζi+1−ζin) n(ζin−ζi−1)
2. 根据上述差分方程，写出求解的fortran 程序（附件1,2）;
3. 根据fortran 程序求出三种格式下的解;
4. 比较各种格式下不同的的数值解，讨论它们的特点。
∆x∆t
四、实验结果分析
利用上述fortran 程序的计算结果，可以观察计算数值，初步判断差分格式的收敛性和稳定性；然后利用EXCEL 或MATLAB 进行绘图，更形象地体现计算结果。以下为利用计算数据在MATLAB 所绘制的图（为了使图更清晰，对比更明显，只画了x ∈[-3,3]的图）：
1.
∆t
∆x
=0.5(
∆x∆t
=2)
A 格式：
B 格式：
C 格式：
从上述三图看出，A 格式收敛不稳定，B 格式不收敛，C 格式收敛稳定，但有数据耗散的情况。总体来说，在
2.
∆t∆x
=0.5的情况下，C 格式计算精度最高。
=1(∆t=1) ∆x
∆t∆x
A 格式：
B 格式：
C 格式：
从上述三图可以看出，A 格式和B 格式依然不收敛，C 格式稳定收敛，没有发生数值耗散，解出的结果最精确。
3.
∆t
∆x
=2(
∆x∆t
=0.5)
A 格式：
B 格式：
C 格式：
从上述三图可以看出，A 、B 、C 格式均不收敛。因为 的值过大，导致时间网格过大，所
∆x∆t
以三种格式都不收敛。
五、实验结论
在为0.5,1时，A 格式收敛不稳定；在为2时，A 不收敛。在0.5,1,2时，B 格式
∆x
∆x
∆x
∆t
∆t
∆t
均不收敛。在而在时，C 格式收敛，但会有数据耗散；在=1时，C 格式稳定收敛；
∆x
∆x
∆t∆t
但在∆x=2时，C 格式不收敛。该结论与收敛性分析和稳定性分析结果一致。
∆t
附件：
（1）Main.f90文件
! 使用说明
!
注意：该程序必须与data.txt 在同一文件夹内
! 若要修改gamma 和delta x，可在data.txt 中分别修改GAMMA 和DX 的值
! 计算完成后请在该程序项目所在文件夹下寻找Output_FormA(B或C).txt 文件查看计算结果 PROGRAM MAIN ! 声明变量
!TN 时间网格数，XN 空间网格数 !DX ，DT/DX(GAMMA) !RB 右边界，LB 左边界 IMPLICIT NONE INTEGER ::TN,XN,I,J REAL ::DX,GAMMA REAL ::LB,RB
REAL , ALLOCATABLE ::U(:,:),X(:)
LOGICAL ALIVE
CHARACTER (len=79)::TEMP CHARACTER (len=2)::FORM ! 声明全局变量
common TN,XN,DX,GAMMA,RB,LB
WRITE (*,*)"*************Computing***********" ! 打开数据文件"data.txt" 录入6个参数 INQUIRE (FILE='data.txt' ,EXIST=ALIVE) IF (ALIVE) THEN
OPEN (UNIT=10,FILE="data.txt" ,STATUS='OLD' ,ACTION='READ' ) READ (10,"(a)")TEMP ! 标题信息
READ (10,*) TN,DX,GAMMA,LB,RB,FORM ! 读取6个参数 WRITE (*,*) "Read successfully" ELSE
WRITE (*,*)"The document doesn't exist." STOP END IF
! 计算空间网格数 XN=(RB-LB)/DX+1
! 确定U 矩阵和X 向量的大小 ALLOCATE (U(XN,TN+1)) ALLOCATE (X(XN)) ! 定义U 初值为0 U=0
! 计算t=0时，U 的函数值U （:,1） DO I=1,XN
X(I)=LB+(I-1)*DX
IF (X(I)>=-1 .AND. X(I)
ELSE IF(X(I)>0 .AND. X(I)
! 调用子程序计算求解
IF (FORM=="A" .or.FORM=="a" ) THEN CALL FormA(U,X)
ELSE IF(FORM=="B" .or.FORM=="b" ) THEN CALL FormB(U,X)
ELSE IF(FORM=="C" .or.FORM=="c" ) THEN CALL FormC(U,X) ELSE
WRITE (*,*)"ERROR:Form fault"
stop
END IF stop END PROGRAM
!***************子程序 A******************** Subroutine FormA(U,X) implicit none INTEGER ::TN,XN,I,J REAL ::DX,GAMMA REAL ::LB,RB
real ::U(XN,TN),X(XN) common TN,XN,DX,GAMMA,RB,LB ! 打开输出文件用于输出数据
OPEN (UNIT=9,FILE="Output_FormA.txt") ! 计算T=TN的值
WRITE (9,*)" 差分格式：格式A" WRITE (9,*)"DT/DX=",GAMMA DO J=2,TN
WRITE (9,*) "TIME=",J-1
WRITE (9,"(A5,10X,A1)") " 坐标X" , "U" DO I=1,XN
IF (I==1)THEN ! 边界点置为相邻点的值 U(I,J)=U(I+1,J-1) ELSE IF(I==XN) then
U(I,J)=U(I-1,J-1) ELSE
U(I,J)=U(I,J-1)-GAMMA*(U(I+1,J-1)-U(I-1,J-1))/2 END IF
WRITE (9,*) X(I),U(I,J) END DO END DO
WRITE (*,*)"***********Finish**********"
write (*,*)" 请在程序文件夹下打开Output_FormA.txt查看结果" return end subroutine
!***************子程序 B******************** Subroutine FormB(U,X) implicit none INTEGER ::TN,XN,I,J REAL ::DX,GAMMA REAL ::LB,RB
real ::U(XN,TN),X(XN) common TN,XN,DX,GAMMA,RB,LB ! 打开输出文件用于输出数据
OPEN (UNIT=9,FILE="Output_FormB.TXT") ! 计算T=TN的值
WRITE (9,*)" 差分格式：格式B" WRITE (9,*)"DT/DX=",GAMMA DO J=2,TN
WRITE (9,*) "TIME=",J-1
WRITE (9,"(A5,10X,A1)") " 坐标X" , "U" DO I=1,XN
IF (I==1)THEN ! 边界点置为相邻点的值 U(I,J)=U(I+1,J-1) ELSE IF(I==XN) then
U(I,J)=U(I-1,J-1) ELSE
U(I,J)=U(I,J-1)-GAMMA*(U(I+1,J-1)-U(I,J-1)) END IF
WRITE (9,"(F6.3,6X,F6.3)") X(I),U(I,J) END DO END DO
WRITE (*,*)"***********Finish**********"
write (*,*)" 请在程序文件夹下打开Output_FormB.txt查看结果" return end subroutine
!***************子程序 C******************** Subroutine FormC(U,X) implicit none INTEGER ::TN,XN,I,J REAL ::DX,GAMMA REAL ::LB,RB
real ::U(XN,TN),X(XN) common TN,XN,DX,GAMMA,RB,LB ! 打开输出文件用于输出数据
OPEN (UNIT=9,FILE="Output_FormC.txt") ! 计算T=TN的值
WRITE (9,*)" 差分格式：格式C" WRITE (9,*)"DT/DX=",GAMMA DO J=2,TN
WRITE (9,*) "TIME=",J-1
WRITE (9,"(A5,10X,A1)") " 坐标X" , "U" DO I=1,XN
IF (I==1)THEN ! 边界点置为相邻点的值
U(I,J)=U(I+1,J-1) ELSE IF(I==XN) then
U(I,J)=U(I-1,J-1) ELSE
U(I,J)=U(I,J-1)-GAMMA*(U(I,J-1)-U(I-1,J-1)) END IF
WRITE (9,"(F6.3,6X,F6.3)") X(I),U(I,J) END DO
write (9,"(80X)") END DO
WRITE (*,*)"***********Finish**********"
write (*,*)" 请在程序文件夹下打开Output_FormC.txt查看结果" return end subroutine
(2)data.txt
TN DX GAMMA LB RB FORM 5 0.05 0.5 -8 8 A
附件说明：主要更改TN,GAMMA 和FORM 的值。分别为时间网格数，DT_DX,差分格式 其中，TN 取任何大于1的整数；GAMMA 可取值为0.5,1,2;FORM 可取为A,B,C
一维对流方程在三种差分格式下的数值解
一、实验目的
用数值方法计算一维对流方程在A 、B 、C 三种差分格式下的解。∆x 取为0.05，∆t为0.5,1,2. 并作相关讨论
ðζðζ+=0, −8≤x≤8 ,t>0 ðtðx
ζ(x,0) =1+x,−1≤x≤0 ζ(x,0) =1−x,0≤x≤1 ζ(x,0) =0,
x1
∆x
二、实验要求
1. 学会对MS-FORTRAN 的基本操作
2. 用Fortran 编写程序计算一维对流方程在A 、B 、C 三种格式下的解。 3. 讨论各种格式下不同的
∆x∆t
三、实验步骤
1. 写出一维对流方程的三种格式的差分方程：
nn
A 格式：ζin+1=ζin+2∆x(ζi+1−ζi−1)
∆t
B 格式：ζin+1=ζin+C 格式：ζin+1=ζin+
∆t∆x∆t∆x
n
(ζi+1−ζin) n(ζin−ζi−1)
2. 根据上述差分方程，写出求解的fortran 程序（附件1,2）;
3. 根据fortran 程序求出三种格式下的解;
4. 比较各种格式下不同的的数值解，讨论它们的特点。
∆x∆t
四、实验结果分析
利用上述fortran 程序的计算结果，可以观察计算数值，初步判断差分格式的收敛性和稳定性；然后利用EXCEL 或MATLAB 进行绘图，更形象地体现计算结果。以下为利用计算数据在MATLAB 所绘制的图（为了使图更清晰，对比更明显，只画了x ∈[-3,3]的图）：
1.
∆t
∆x
=0.5(
∆x∆t
=2)
A 格式：
B 格式：
C 格式：
从上述三图看出，A 格式收敛不稳定，B 格式不收敛，C 格式收敛稳定，但有数据耗散的情况。总体来说，在
2.
∆t∆x
=0.5的情况下，C 格式计算精度最高。
=1(∆t=1) ∆x
∆t∆x
A 格式：
B 格式：
C 格式：
从上述三图可以看出，A 格式和B 格式依然不收敛，C 格式稳定收敛，没有发生数值耗散，解出的结果最精确。
3.
∆t
∆x
=2(
∆x∆t
=0.5)
A 格式：
B 格式：
C 格式：
从上述三图可以看出，A 、B 、C 格式均不收敛。因为 的值过大，导致时间网格过大，所
∆x∆t
以三种格式都不收敛。
五、实验结论
在为0.5,1时，A 格式收敛不稳定；在为2时，A 不收敛。在0.5,1,2时，B 格式
∆x
∆x
∆x
∆t
∆t
∆t
均不收敛。在而在时，C 格式收敛，但会有数据耗散；在=1时，C 格式稳定收敛；
∆x
∆x
∆t∆t
但在∆x=2时，C 格式不收敛。该结论与收敛性分析和稳定性分析结果一致。
∆t
附件：
（1）Main.f90文件
! 使用说明
!
注意：该程序必须与data.txt 在同一文件夹内
! 若要修改gamma 和delta x，可在data.txt 中分别修改GAMMA 和DX 的值
! 计算完成后请在该程序项目所在文件夹下寻找Output_FormA(B或C).txt 文件查看计算结果 PROGRAM MAIN ! 声明变量
!TN 时间网格数，XN 空间网格数 !DX ，DT/DX(GAMMA) !RB 右边界，LB 左边界 IMPLICIT NONE INTEGER ::TN,XN,I,J REAL ::DX,GAMMA REAL ::LB,RB
REAL , ALLOCATABLE ::U(:,:),X(:)
LOGICAL ALIVE
CHARACTER (len=79)::TEMP CHARACTER (len=2)::FORM ! 声明全局变量
common TN,XN,DX,GAMMA,RB,LB
WRITE (*,*)"*************Computing***********" ! 打开数据文件"data.txt" 录入6个参数 INQUIRE (FILE='data.txt' ,EXIST=ALIVE) IF (ALIVE) THEN
OPEN (UNIT=10,FILE="data.txt" ,STATUS='OLD' ,ACTION='READ' ) READ (10,"(a)")TEMP ! 标题信息
READ (10,*) TN,DX,GAMMA,LB,RB,FORM ! 读取6个参数 WRITE (*,*) "Read successfully" ELSE
WRITE (*,*)"The document doesn't exist." STOP END IF
! 计算空间网格数 XN=(RB-LB)/DX+1
! 确定U 矩阵和X 向量的大小 ALLOCATE (U(XN,TN+1)) ALLOCATE (X(XN)) ! 定义U 初值为0 U=0
! 计算t=0时，U 的函数值U （:,1） DO I=1,XN
X(I)=LB+(I-1)*DX
IF (X(I)>=-1 .AND. X(I)
ELSE IF(X(I)>0 .AND. X(I)
! 调用子程序计算求解
IF (FORM=="A" .or.FORM=="a" ) THEN CALL FormA(U,X)
ELSE IF(FORM=="B" .or.FORM=="b" ) THEN CALL FormB(U,X)
ELSE IF(FORM=="C" .or.FORM=="c" ) THEN CALL FormC(U,X) ELSE
WRITE (*,*)"ERROR:Form fault"
stop
END IF stop END PROGRAM
!***************子程序 A******************** Subroutine FormA(U,X) implicit none INTEGER ::TN,XN,I,J REAL ::DX,GAMMA REAL ::LB,RB
real ::U(XN,TN),X(XN) common TN,XN,DX,GAMMA,RB,LB ! 打开输出文件用于输出数据
OPEN (UNIT=9,FILE="Output_FormA.txt") ! 计算T=TN的值
WRITE (9,*)" 差分格式：格式A" WRITE (9,*)"DT/DX=",GAMMA DO J=2,TN
WRITE (9,*) "TIME=",J-1
WRITE (9,"(A5,10X,A1)") " 坐标X" , "U" DO I=1,XN
IF (I==1)THEN ! 边界点置为相邻点的值 U(I,J)=U(I+1,J-1) ELSE IF(I==XN) then
U(I,J)=U(I-1,J-1) ELSE
U(I,J)=U(I,J-1)-GAMMA*(U(I+1,J-1)-U(I-1,J-1))/2 END IF
WRITE (9,*) X(I),U(I,J) END DO END DO
WRITE (*,*)"***********Finish**********"
write (*,*)" 请在程序文件夹下打开Output_FormA.txt查看结果" return end subroutine
!***************子程序 B******************** Subroutine FormB(U,X) implicit none INTEGER ::TN,XN,I,J REAL ::DX,GAMMA REAL ::LB,RB
real ::U(XN,TN),X(XN) common TN,XN,DX,GAMMA,RB,LB ! 打开输出文件用于输出数据
OPEN (UNIT=9,FILE="Output_FormB.TXT") ! 计算T=TN的值
WRITE (9,*)" 差分格式：格式B" WRITE (9,*)"DT/DX=",GAMMA DO J=2,TN
WRITE (9,*) "TIME=",J-1
WRITE (9,"(A5,10X,A1)") " 坐标X" , "U" DO I=1,XN
IF (I==1)THEN ! 边界点置为相邻点的值 U(I,J)=U(I+1,J-1) ELSE IF(I==XN) then
U(I,J)=U(I-1,J-1) ELSE
U(I,J)=U(I,J-1)-GAMMA*(U(I+1,J-1)-U(I,J-1)) END IF
WRITE (9,"(F6.3,6X,F6.3)") X(I),U(I,J) END DO END DO
WRITE (*,*)"***********Finish**********"
write (*,*)" 请在程序文件夹下打开Output_FormB.txt查看结果" return end subroutine
!***************子程序 C******************** Subroutine FormC(U,X) implicit none INTEGER ::TN,XN,I,J REAL ::DX,GAMMA REAL ::LB,RB
real ::U(XN,TN),X(XN) common TN,XN,DX,GAMMA,RB,LB ! 打开输出文件用于输出数据
OPEN (UNIT=9,FILE="Output_FormC.txt") ! 计算T=TN的值
WRITE (9,*)" 差分格式：格式C" WRITE (9,*)"DT/DX=",GAMMA DO J=2,TN
WRITE (9,*) "TIME=",J-1
WRITE (9,"(A5,10X,A1)") " 坐标X" , "U" DO I=1,XN
IF (I==1)THEN ! 边界点置为相邻点的值
U(I,J)=U(I+1,J-1) ELSE IF(I==XN) then
U(I,J)=U(I-1,J-1) ELSE
U(I,J)=U(I,J-1)-GAMMA*(U(I,J-1)-U(I-1,J-1)) END IF
WRITE (9,"(F6.3,6X,F6.3)") X(I),U(I,J) END DO
write (9,"(80X)") END DO
WRITE (*,*)"***********Finish**********"
write (*,*)" 请在程序文件夹下打开Output_FormC.txt查看结果" return end subroutine
(2)data.txt
TN DX GAMMA LB RB FORM 5 0.05 0.5 -8 8 A
附件说明：主要更改TN,GAMMA 和FORM 的值。分别为时间网格数，DT_DX,差分格式 其中，TN 取任何大于1的整数；GAMMA 可取值为0.5,1,2;FORM 可取为A,B,C