扩散方程的差分解法

5.3计算程序
5.3.1显格式
!----------------------------------------显格式求解扩散方程---------------------------------------------------!
!----------------------参数设置-----------------------------------------！
!---------------------初始条件-----------------------------------------!
do j=1,nx
x=(j-1)*dx
if (x.le.sx/2) then
u(j,1)=2*x /sx
else
u(j,1)=2*(1-x/sx)
endif
enddo
2.有限差分法
有限差分法是数值计算解微分方程古老的方法之一，也是系统化地、数值地求解数学物理方法的方程。其控制方程中的导数用离散点上函数值的差商代替。
差分格式可以分为显格式和隐格式。所谓显格式是指在任一结点上因变量在新是时间层上的值可以通过之前的时间层上相邻结点变量的值显式解出来。由于这些层的变量值是已知的，当时间向前推进时，空间点上的新的变量值就只需逐点计算就行了，因此显格式计算起来比较省事。隐格式则是指任一结点上变量在新的时间层的值，不能通过之前的时间层上相邻结点的值显式解出来，它不仅与之前的时间层上的已知值有关，而且也与新时间层的相邻结点的变量值有关。因而一个差分方程常常包括几个相邻结点上的未知数，未知数的个数取决于格式的构成形式。为了解出这些未知数需要联立新的方程，而每引进一个新的方程往往又同时引进了新的未知数。因此，隐格式总是伴随着求解巨大的代数方程组。隐格式的主要缺点是计算工作量大，因而不如显格式计算得快，但这只是就时间步长一样的情况而言的。隐格式的主要优点是时间步长可以比显格式能够采用的最大步长大很多。显格式的时间步长受到稳定性条件的限制，而隐格式则几乎不受限制。
对一维扩散方程：
（7）
采用显格式差分格式，令 ， ，则
（8）
用TaLeabharlann Baidulor级数展开代入上述差分方程中，则有
（9）
当 ， 时，上式化为
（10）
可见，此差分格式所构成的差分方程与原来的微分方程是相容的，故该显格式为相容格式。
采用全隐格式差分格式，令 ， ，则
（11）
用Taylor级数展开代入上述差分方程中，则有
enddo
!-----------!
!-------------------每隔15s输出沿程变化------------------------!
do n=1,nt
t=(n-1)*dt
if (mod(int(t*1000),15000).eq.0) then
write(1,*) 't=',t,'s'
a(1)=1
do i=1,nx-1
d(i)=-r
c(i)=-r
a(i+1)=1+2*r
enddo
c(1)=0
!-----------!
!-----------------------------三对角方程组常数项-------------------------!
do i=1,nx
b(i)=u(i,1)
扩散方程的差分解法
在研究热传导过程、扩散过程、边界层现象时，我们常常遇到抛物型方程，这类方程中最典型、最简单的就是热传导方程。热传导方程中的自变量中包括时间t，它是描述一种随时间变化的物理过程，即所谓不定常现象。这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。
do j=1,nx
write(1,*) (j-1)*dx,u(j,n)
enddo
endif
enddo
!-----------!
!-------------------每隔2.5m输出随时间变化------------------------!
do j=1,nx
x=dx*(j-1)
if (mod(int(x*10),25).eq.0) then
3.方程的离散
3.1显格式
采用时间前差及第n时间层的空间中心差，得一维扩散方程的显格式解：
（4）
即
（5）
式中，
3.2全隐格式
采用时间前差及第（n+1）时间层的空间中心差，得一维扩散方程的全隐格式解：
（6）
4.差分解的基本问题
差分解的基本问题包括：适定性、相容性、收敛性和稳定性四个方面。
4.1适定性
在用差分方程作微分方程数值解时，首先，要求微分方程的问题是适定的。所谓适定性问题是指这一微分方程在一定的初始和边界条件下要有唯一解，并且在初始条件和边界条件稍有改变时，微分方程的解也只是稍有偏离。从数学的角度讲，若微分方程的解存在并且是唯一的，同时连续依赖于数据（初始条件、边界条件），则问题是适定的。
（12）
当 ， 时，上式化为
（13）
可见，此差分格式所构成的差分方程与原来的微分方程是相容的，故该全隐格式为相容格式。
4.3稳定性
差分法计算中所产生的误差（舍入误差，参数误差等）随时间衰减或不增大，则称离散格式是稳定的，反之，则是不稳定的。分析差分方程稳定性有不同的方法，如矩阵方法，谐波分析法等。下面用谐波分析法对以上两种离散格式的稳定性进行分析：
do j=1,nx
x=(j-1)*dx
if (x.le.sx/2) then
u(j,1)=2*x/sx
else
u(j,1)=2*(1-x/sx)
endif
enddo
!-----------!
!-----------------------------三对角方程组系数矩阵-------------------------!
r=alfa*dt/dx**2
!----------------------边界条件----------------------------------------!
do n=1,nt
u(1,n)=0
u(nx,n)=0
enddo
!----------!
!---------------------初始条件-----------------------------------------!
由以上对一维扩散问题的分析，可知，求解一维扩散方程需给定初始条件及边界条件。
在本文计算中，取 ， 。
初始条件（ 时）
（29）
边界条件为
（30）
其初始时刻（ ）时的u分布如图1所示，x=0m处u随时间变化情况如图2所示，x=10m处u随时间变化情况如图3所示。
图1初始时刻u分布图
图2 x=0处u随时间变化图
!----------------------参数设置-----------------------------------------！
!implicit double precision (a-h,o-z)
parameter (nt=20001,nx=101) !nt:时间节点数;nx:空间节点数
（28）
则称差分格式是收敛的。
拉克斯（Lax）等价定理指出，如果问题是适定的，并且差分格式满足相容性条件，那么差分格式的稳定性就是该格式收敛性的充分而必要的条件。
由4.1、4.2和4.3的分析可知，显格式和全隐格式都满足适定性、相容性及稳定性的条件，因而这两种格式满足收敛性要求。
5.差分方程的求解
5.1初始条件及边界条件
在本文的一维扩散方程中，给定了初始条件，同时在区域的左、右端边界给定了边界条件，满足适定性要求。
4.2相容性
将一个偏微分方程用差分格式化为相应差分方程，当步长 和 趋近于零时，这个差分方程应当收敛于原微分方程，也就是说，相应的差分方程和微分方程之间的截断误差在任一时刻任一网格点上均应趋近于零，这样的差分方程和微分方程才是相容的。
本文主要运用有限差分法对一维扩散方程进行求解，并对差分解的适定性、相容性、收敛性及稳定性进行分析，同时与解析解进行对比。
1.扩散方程
一维扩散方程为：
（1）
式中，u为因知量， 为扩散系数，x为坐标，t为时间。
其定解条件如下：
初始条件：
（2）
边界条件：
（3）
一般假定函数 ， ， 满足连接条件，即 ， 。
图3 x=10m处u随时间变化图
5.2空间步长及时间步长
通过对差分格式的稳定性分析知，显格式空间步长与时间步长间应满足一定的关系，即式（21）。本文选用的步长如表1所示，分别用显格式和全隐格式进行数值计算，差分网格如图4所示。
表1时间步长及空间步长取值表
时间步长 /s
空间步长 /m
时间步数
计算时长/s
write(2,*) 'x=',x,'m'
do n=1,nt,200
write(2,*) (n-1)*dt,u(j,n)
enddo
endif
enddo
!-----------!
end
5.3.2全隐格式
!----------------------------------------全隐格式求解扩散方程-----------------------------------------------!
（1）显格式
显格式的差分方程为
（14）
即
（15）
其误差方程为
（16）
任取一k次谐波分量
（17）
则
（18）
误差放大因子为
（19）
要满足稳定性条件，则要求对所有的k值均有 ，须 ，即 。
因此，一维扩散方程显格式的稳定性条件为
（20）
（2）全隐格式
全隐格式的差分方程为
（21）
即
（22）
其误差方程为
（23）
对于抛物线方程，这是初值问题，要求给出初始条件，并且在区域边界上给出边界条件。在本文的一维扩散问题中，除了给出t=0时的初值 外，应在左、右两端边界，即 及 ，各给定一个边界条件，第一类边界条件u的值，或给定第二类边界条件 的值或给定第三类边界条件u与 的组合。
由于抛物型方程具有扩散性质，它往往是随时间扩散的，u值将不断因扩散而衰减。
enddo
!-----------!
!----------------------------------求解差分方程-----------------------------!
do n=2,nt
call systri(nx,d,a,c,b,xx)
parameter (nt=20001,nx=101)!nt:时间节点数;nx:空间节点数
dimension u(nx,nt)
Sx=10.0!计算长度（m）
dt=0.003!时间步长（s）
dx=0.1!空间步长（m）
alfa=1.0!扩散系数
!-----------!
open(1,file='显格式沿程变化.txt')
dimension u(nx,nt),a(nx),b(nx),c(nx-1),d(nx-1),xx(nx)
Sx=10!计算长度（m）
dt=0.003!时间步长（s）
dx=0.1!空间步长（m）
alfa=1.0!扩散系数
!-----------!
open(1,file='全隐格式沿程变化.txt')
open(2,file='全隐格式随时间变化.txt')
空间步数
0.003
0.1
500000
1500
100
0.3
图4差分网格示意图
5.2求解方法
对于显格式，有 ，若n时间层上的u已知，则可直接求解出（n+1）时间层上的u值。因此，由于给出了初始条件，显格式无需迭代，直接从t=0时刻开始，逐一计算下一个时间层上的u值，便可求解出各时间层上各空间点的u值。
对于全隐格式，由于 与 ， 有联系，不能直接求解，必须联解代数方程组。此时方程组为三对角方程组，可采用追赶法进行求解。
open(2,file='显格式随时间变化.txt')
r=alfa*dt/dx**2
!----------------------边界条件----------------------------------------!
do n=1,nt
u(1,n)=0
u(nx,n)=0
enddo
!----------!
任取一k次谐波分量
（24）
则
， （25）
则误差方程为
（26）
误差放大因子为
（27）
要满足稳定性条件，则要求对所有的k值均有 。从（28）式中可以看出，当 （即 ）时， 恒成立。因此，全隐格式是无条件稳定的。
4.4收敛性
如果差分方程的解为 ，微分方程的解为 ，若当 ， 时，差分方程的解与微分方程的解之差
!-----------!
!--------------------求解差分方程------------------------------------!
do n=2,nt
do j=2,nx-1
u(j,n)=u(j,n-1)+r*(u(j-1,n-1)-2*u(j,n-1)+u(j+1,n-1))
enddo