扩散方程的差分解法
扩散方程的数值差分解法
扩散方程的数值差分解法作者：刘浩庭来源：《价值工程》2019年第29期摘要：扩散现象是其初始密度不均匀分布引起的，它会对物质粒子的分布状态产生影响，最终达到物质在空间均匀分布状态。
本研究通过分离变量法对给定条件的粒子浓度在一维空间分布下的扩散现象的解析解进行计算，同时结合使用欧拉法利用计算物理的方法对扩散过程进行了数值模拟。
通过对比理论数据与模拟实验数据对等离子体一维扩散现象进行阐释与讨论。
Abstract： The diffusion phenomenon is caused by the difference of the initial density. It affects the motion of the particles and finally make all the particles into the uniformly distribution. In this study， the analytical solution of the diffusion phenomenon of a given particle concentration in a one-dimensional space is calculated by using the method of separation of variables， and the diffusion process is simulated numerically by using the computational physics method combined with Euler's polygonal arc method. The one-dimensional plasma diffusion phenomenon is explained by comparing the theoretical data with the simulated experimental data.關键词：一维扩散;分离变量法;欧拉法Key words： one-dimensional diffusion;the method of separation of variables;Euler's polygonal arc method中图分类号：O122.2; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;文献标识码：A; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; 文章编号：1006-4311（2019）29-0272-041; 简介受控热核聚变是受世人瞩目的前沿重要课题，其目的便是探索清洁可持续的新能源。
求解空间分数阶扩散方程和对流扩散方程的有限差分格式研究
求解空间分数阶扩散方程和对流扩散方程的有限差分格式研究1 空间分数阶扩散方程有限差分格式研究空间分数阶扩散方程是一类非线性偏微分方程，广泛应用于化学、生物、地理、物理等领域的模拟和研究中。
由于其阶数为分数阶，因此其求解方法与常规的整数阶偏微分方程有所不同。
##1.1 基本方程及边值条件空间分数阶扩散方程基本形式为：$$\frac{\partial^{\alpha}u}{\partial t^{\alpha}}=D\frac{\partial^2u}{\partial x^2}$$ 其中，$0<\alpha<1$为分数阶，$D$为扩散系数，$u(x,t)$为扩散物体在空间$x$和时间$t$的浓度分布。
边值条件通常为：$$u(x,0)=f(x)$$$$u(0,t)=u(L,t)=0$$其中，$f(x)$为初始浓度分布，$L$为空间长度。
##1.2 有限差分格式为了在计算机上求解空间分数阶扩散方程，需要将其离散化为有限差分格式。
常用的有限差分格式为Caputo分数阶导数格式和Grünwald-Letnikov分数阶导数格式。
这里以Caputo分数阶导数格式为例，其形式为：$$\frac{\partial^{\alpha}u}{\partial t^{\alpha}}\approx\frac{1}{\Gamma(1-\alpha)}\int_0^t\frac{\partial u}{\partial s}(t-s)^{-\alpha}ds$$$$\frac{\partial^2u}{\partial x^2}\approx\frac{u(x+\Delta x)-2u(x)+u(x-\Deltax)}{\Delta x^2}$$将上述两式带入空间分数阶扩散方程中，得到：$$\frac{1}{\Gamma(1-\alpha)}\int_0^t\frac{\partial u(x,s)}{\partial s}(t-s)^{-\alpha}ds=D\frac{u(x+\Delta x)-2u(x)+u(x-\Delta x)}{\Delta x^2}$$可得到迭代公式：$$u_i^{n+1}=\frac{(1-\theta)\Delta t^{\alpha}}{\Gamma(2-\alpha)\Deltax^2}u_{i+1}^n+\frac{2\theta\Delta t^{\alpha}}{\Gamma(2-\alpha)\Deltax^2}u_i^{n+1/2}+\frac{(1-\theta)\Delta t^{\alpha}}{\Gamma(2-\alpha)\Delta x^2}u_{i-1}^n+\frac{\Delta t}{\Delta x^2}f_i^n$$其中，$u_i^n$表示在$x=i\Delta x$、$t=n\Delta t$时的浓度值，$f_i^n$表示边界条件。
扩散方程的差分解法
扩散方程的差分解法在研究热传导过程、扩散过程、边界层现象时，我们常常遇到抛物型方程，这类方程中最典型、最简单的就是热传导方程。
热传导方程中的自变量中包括时间t ，它是描述一种随时间变化的物理过程，即所谓不定常现象。
这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。
本文主要运用有限差分法对一维扩散方程进行求解，并对差分解的适定性、相容性、收敛性及稳定性进行分析，同时与解析解进行对比。
1.扩散方程一维扩散方程为：22u u t xα∂∂=∂∂ （1）式中，u 为因知量，α为扩散系数，x 为坐标，t 为时间。
其定解条件如下： 初始条件： (,0)() 0x u x f x L =≤≤（2）边界条件： 12(0,)() , (,)()u t f t u L t f t ==（3） 一般假定函数()f x ，1()f t ，2()f t 满足连接条件，即1(0)(0) f f =，2()(0) f L f =。
2.有限差分法有限差分法是数值计算解微分方程古老的方法之一，也是系统化地、数值地求解数学物理方法的方程。
其控制方程中的导数用离散点上函数值的差商代替。
差分格式可以分为显格式和隐格式。
所谓显格式是指在任一结点上因变量在新是时间层上的值可以通过之前的时间层上相邻结点变量的值显式解出来。
由于这些层的变量值是已知的，当时间向前推进时，空间点上的新的变量值就只需逐点计算就行了，因此显格式计算起来比较省事。
隐格式则是指任一结点上变量在新的时间层的值，不能通过之前的时间层上相邻结点的值显式解出来，它不仅与之前的时间层上的已知值有关，而且也与新时间层的相邻结点的变量值有关。
因而一个差分方程常常包括几个相邻结点上的未知数，未知数的个数取决于格式的构成形式。
为了解出这些未知数需要联立新的方程，而每引进一个新的方程往往又同时引进了新的未知数。
因此，隐格式总是伴随着求解巨大的代数方程组。
一维非稳态对流-扩散方程的隐式中心差分格式
一维非稳态对流-扩散方程的隐式中心差分格式非稳态对流-扩散方程：$$\frac{\partial u}{\partial t}+\frac{\partial}{\partialx}\left(F_{c}\left(u\right)\right)=\underbrace{\frac{\partial}{\ partial x}\left(V_{x}\left(x\right)\frac{\partial u}{\partial x}\right)}_{扩散项}+S\left(x,t\right)$$其中，$u\left(x,t\right)$为温度、浓度等物理量，$F_{c}\left(u\right)$为流动的拖曳力，$V_{x}\left(x\right)$为扩散系数，$S\left(x,t\right)$为源项。
考虑前向差分格式，$i,j=\cdots,-2,-1,0,1,2,\cdots$表示位置，$n$表示时刻：$$u_{i}^{n+1}=u_{i}^{n}-\frac{\Delta t}{\Deltax}\left[F_{c}\left(u_{i+1}^{n}\right)-F_{c}\left(u_{i}^{n}\right)\right]+\frac{\Delta t}{2\Deltax}\left[V_{i+1/2}^{n+1/2}\left(u_{i+1}^{n}-u_{i}^{n}\right)+V_{i-1/2}^{n+1/2}\left(u_{i}^{n}-u_{i-1}^{n}\right)\right]+\Delta t\cdot S_{i}^{n}$$其中，$V_{i+1/2}^{n+1/2}=V_{x}\left(x_{i+1/2}^{n+1/2}\right)$，$x_{i+1/2}^{n+1/2}=\frac{1}{2}\left(x_{i+1}^{n}+x_{i+1}^{n+1}\ri ght)$。
diffusion扩散模型运用的算法
diffusion扩散模型运用的算法Diffusion扩散模型是一种用于研究物质在空间中传播和扩散的数学模型。
它可以描述分子、热量、能量等在不同浓度或温度下的自然扩散现象。
该模型广泛应用于物理学、化学、生物学等领域，并被用于解决各种实际问题。
在扩散模型中，物质的传播可以通过扩散方程来描述。
扩散方程是一个偏微分方程，它描述了物质在空间中的浓度随时间的变化。
该方程的形式如下：∂C/∂t = D∇²C其中，C是物质的浓度，t是时间，D是扩散系数，∇²是拉普拉斯算子。
这个方程表明，物质浓度的变化率等于扩散系数乘以浓度的二阶空间导数。
在实际应用中，为了解决扩散方程，可以采用不同的算法。
下面介绍两种常用的算法：有限差分法和有限元法。
有限差分法是一种将连续方程离散化为差分方程的方法。
它将空间和时间分成若干个小区间，然后用差分近似代替微分，从而将连续方程转化为离散方程。
在扩散模型中，可以将空间划分为网格点，然后根据差分近似计算每个网格点的浓度。
通过迭代计算，可以得到整个空间中物质浓度的分布。
有限元法是一种将连续方程离散化为有限个元素方程的方法。
它将空间划分为若干个小单元，然后用一组基函数逼近每个小单元内的物质浓度。
通过求解元素方程，可以得到整个空间中物质浓度的近似解。
有限元法相对于有限差分法具有更高的精度和灵活性，适用于复杂的几何形状和边界条件。
除了有限差分法和有限元法，还有其他一些算法可以用于解决扩散模型。
例如，蒙特卡洛方法可以通过随机模拟分子运动来模拟扩散过程。
这种方法基于概率思想，通过大量的模拟实验来估计物质浓度的分布。
蒙特卡洛方法不依赖于方程的解析解，适用于复杂的非线性和非均匀问题。
扩散模型的算法应用范围广泛。
在物理学中，可以用扩散模型来研究热传导、电子输运等现象。
在化学中，可以用扩散模型来研究溶质在溶液中的传输和反应。
在生物学中，可以用扩散模型来研究细胞内物质的传输和扩散。
Diffusion扩散模型是一种重要的数学模型，通过不同的算法可以解决各种实际问题。
有限差分法求解扩散方程的步骤
有限差分法求解扩散方程的步骤有限差分法是求解扩散方程的一种有效方法，简称FDM，有限差分法能够解决复杂的扩散方程，可以看作数值计算在扩散方程中的一个应用。
一般情况下，有限差分法求解扩散方程是通过将扩散方程分解为两部分：非线性问题和线性问题，分别用不同的求解方法解决。
在这篇文章中，我们将讨论使用有限差分法求解扩散方程的步骤，帮助读者更好地理解有限差分法。
第一步：建立数值解模型。
有限差分法求解扩散方程，首先要建立数值解模型。
可以将扩散方程的区域划分为若干个小矩形，用每个小矩形的中心的值代表这一区域的大致状态，然后计算每一部分的有限差分，从而建立起数值解模型。
第二步：求解线性问题。
这一步用来求解扩散方程中的线性部分，包括：首先，对离散点的值进行定义；其次，在离散点之间建立差分关系；最后，根据上述关系，求解离散点的值。
第三步：求解非线性问题。
有限差分法还可以求解扩散方程中的非线性部分。
可以先将非线性部分转化为线性部分，然后求解，也可以使用迭代法求解。
第四步：检查模型的正确性。
有限差分法求解扩散方程后，需要检查模型的正确性，可以使用数值积分方法、定性方法或定量分析等方法来检查求解结果的正确性。
总之，有限差分法求解扩散方程的五个步骤是：建立数值解模型，求解线性问题，求解非线性问题，进行模型校正以及检查模型的正确性。
在这五个步骤中，第一步特别重要，因为它是整个有限差分法求解过程的基础，如果第一步建立的模型不合理，就不可能得到准确的结果。
有限差分法的运用不仅当前广泛，而且在未来也有很大的发展前景。
由于有限差分法求解扩散方程的步骤具有一定的复杂性，因此有必要在深入研究有限差分法求解扩散方程之前，充分理解这一步骤。
综上所述，有限差分法求解扩散方程的步骤是：建立数值解模型，求解线性问题，求解非线性问题，进行模型校正以及检查模型的正确性。
有限差分法在求解扩散问题方面具有一定的优势，适用范围也较广，此外还有很大的发展前景。
一类二维稳态对流——扩散方程的有限差分法
一类二维稳态对流——扩散方程的有限差分法一维稳态扩散方程描述了物质在一维空间中的扩散行为。
然而，在某些情况下，我们需要研究物质在二维平面中的扩散行为，例如热传导、流体传输等。
本文将介绍一类二维稳态对流-扩散方程的有限差分法。
二维稳态对流-扩散方程可以写作：∇·(D∇u) + ∇·(cu) + fu = 0 —— (1)其中，D是扩散系数，c是速度场，u是待求解的物理量，f是源项。
在这个方程中，第一项表示物质的扩散项，第二项表示对流项，第三项表示源项。
我们需要求解方程(1)，找到u的分布。
为了应用有限差分法来求解二维稳态对流-扩散方程，需要将二维空间离散化为一个网格。
假设我们将x方向离散为Nx个等距的节点，y方向离散为Ny个等距的节点，那么我们可以得到一个(Nx+1)×(Ny+1)的网格。
我们在网格节点上定义未知量u，然后将方程(1)对节点处的u进行离散化。
首先，我们对方程(1)的扩散项进行离散化。
我们使用五点差分格式来近似二维Laplace算符∇·(D∇u)。
对于网格节点(x,y)，我们可以得到以下差分格式：(Dij(xi+1,yj)ui+1,j + Dij(xi-1,yj)ui-1,j +Dij(xi,yj+1)ui,j+1 + Dij(xi,yj-1)ui,j-1 -4Dij(xi,yj)ui,j) / ∆x^2 + (Dij(xi,yj)ui,j) / ∆y^2其中，∆x和∆y是网格步长，Dij是扩散系数。
接下来，我们对方程(1)的对流项进行离散化。
我们使用中心差分格式来近似二维梯度算符∇·(cu)。
对于网格节点(x,y)，我们可以得到以下差分格式：(cxi+1/2,yj(ui+1,j - ui,j)) / ∆x + (cxi-1/2,yj(ui,j - ui-1,j)) / ∆x + (cyi,j+1/2(ui,j+1 - ui,j)) / ∆y + (cyi,j-1/2(ui,j - ui,j-1)) / ∆y其中，cxi+1/2,yj、cxi-1/2,yj、cyi,j+1/2和cyi,j-1/2是速度场在节点(x,y)处的中心点处的x和y分量。
时间分数阶扩散方程的三次样条差分格式
时间分数阶扩散方程的三次样条差分格式时间分数阶扩散方程的三次样条差分格式如下：$$。
\frac{\partial^{\alpha}u(x_i, t^n)}{\partial t^{\alpha}} = D \frac{\partial^2u(x_i, t^n)}{\partial x^2}, \quad 0 < \alpha < 1。
$$。
其中，$\alpha$为时间分数阶，$D$为扩散系数，$u(x_i, t^n)$为在$x_i$处和$t^n$时刻的解。
为方便表示，下面将$t^n$简单表示为$t$。
采用前向差分和后向差分的一阶导数定义，可以得到时间导数的差分格式：$$。
\frac{\partial^{\alpha}u(x_i, t)}{\partial t^{\alpha}}\approx \frac{u(x_i, t) - u(x_i, t - \Delta t)}{\Deltat^{\alpha}}, \quad 0 < \alpha < 1。
$$。
其中，$\Delta t$为时间步长。
由于采用了前向差分，此格式为显式差分格式。
若采用三点中心差分法近似求解二阶导数，则。
$$。
\frac{\partial^2u(x_i, t)}{\partial x^2} \approx\frac{u(x_{i+1},t) - 2u(x_i,t) + u(x_{i-1},t)}{\Delta x^2}。
$$。
其中，$\Delta x$为空间步长。
将上式代入原方程得到：$$。
\frac{u(x_i, t) - u(x_i, t - \Delta t)}{\Delta t^{\alpha}} = D \frac{u(x_{i+1},t) - 2u(x_i,t) + u(x_{i-1},t)}{\Delta x^2}。
$$。
再将时间步长和空间步长合并，得到：$$。
u_i^n - u_i^{n-1} = D\Delta t^{\alpha} \frac{u_{i+1}^n -2u_i^n + u_{i-1}^n}{\Delta x^2}。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
（24）
则
， （25）
则误差方程为
（26）
误差放大因子为
（27）
要满足稳定性条件，则要求对所有的k值均有 。从（28）式中可以看出，当 （即 ）时， 恒成立。因此，全隐格式是无条件稳定的。
4.4收敛性
如果差分方程的解为 ，微分方程的解为 ，若当 ， 时，差分方程的解与微分方程的解之差
扩散方程的差分解法
在研究热传导过程、扩散过程、边界层现象时，我们常常遇到抛物型方程，这类方程中最典型、最简单的就是热传导方程。热传导方程中的自变量中包括时间t，它是描述一种随时间变化的物理过程，即所谓不定常现象。这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。
write(2,*) 'x=',x,'m'
do n=1,nt,200
write(2,*) (n-1)*dt,u(j,n)
enddo
if
enddo
!-----------!
end
5.3.2全隐格式
!----------------------------------------全隐格式求解扩散方程-----------------------------------------------!
由以上对一维扩散问题的分析，可知，求解一维扩散方程需给定初始条件及边界条件。
在本文计算中，取 ， 。
初始条件（ 时）
（29）
边界条件为
（30）
其初始时刻（ ）时的u分布如图1所示，x=0m处u随时间变化情况如图2所示，x=10m处u随时间变化情况如图3所示。
图1初始时刻u分布图
图2 x=0处u随时间变化图
!---------------------初始条件-----------------------------------------!
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
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')
2.有限差分法
有限差分法是数值计算解微分方程古老的方法之一，也是系统化地、数值地求解数学物理方法的方程。其控制方程中的导数用离散点上函数值的差商代替。
差分格式可以分为显格式和隐格式。所谓显格式是指在任一结点上因变量在新是时间层上的值可以通过之前的时间层上相邻结点变量的值显式解出来。由于这些层的变量值是已知的，当时间向前推进时，空间点上的新的变量值就只需逐点计算就行了，因此显格式计算起来比较省事。隐格式则是指任一结点上变量在新的时间层的值，不能通过之前的时间层上相邻结点的值显式解出来，它不仅与之前的时间层上的已知值有关，而且也与新时间层的相邻结点的变量值有关。因而一个差分方程常常包括几个相邻结点上的未知数，未知数的个数取决于格式的构成形式。为了解出这些未知数需要联立新的方程，而每引进一个新的方程往往又同时引进了新的未知数。因此，隐格式总是伴随着求解巨大的代数方程组。隐格式的主要缺点是计算工作量大，因而不如显格式计算得快，但这只是就时间步长一样的情况而言的。隐格式的主要优点是时间步长可以比显格式能够采用的最大步长大很多。显格式的时间步长受到稳定性条件的限制，而隐格式则几乎不受限制。
!-----------!
!--------------------求解差分方程------------------------------------!
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
在本文的一维扩散方程中，给定了初始条件，同时在区域的左、右端边界给定了边界条件，满足适定性要求。
4.2相容性
将一个偏微分方程用差分格式化为相应差分方程，当步长 和 趋近于零时，这个差分方程应当收敛于原微分方程，也就是说，相应的差分方程和微分方程之间的截断误差在任一时刻任一网格点上均应趋近于零，这样的差分方程和微分方程才是相容的。
3.方程的离散
3.1显格式
采用时间前差及第n时间层的空间中心差，得一维扩散方程的显格式解：
（4）
即
（5）
式中，
3.2全隐格式
采用时间前差及第（n+1）时间层的空间中心差，得一维扩散方程的全隐格式解：
（6）
4.差分解的基本问题
差分解的基本问题包括：适定性、相容性、收敛性和稳定性四个方面。
4.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
!-----------!
!-----------------------------三对角方程组系数矩阵-------------------------!
（1）显格式
显格式的差分方程为
（14）
即
（15）
其误差方程为
（16）
任取一k次谐波分量
（17）
则
（18）
误差放大因子为
（19）
要满足稳定性条件，则要求对所有的k值均有 ，须 ，即 。
因此，一维扩散方程显格式的稳定性条件为
（20）
（2）全隐格式
全隐格式的差分方程为
（21）
即
（22）
其误差方程为
（23）
本文主要运用有限差分法对一维扩散方程进行求解，并对差分解的适定性、相容性、收敛性及稳定性进行分析，同时与解析解进行对比。
1.扩散方程
一维扩散方程为：
（1）
式中，u为因知量， 为扩散系数，x为坐标，t为时间。
其定解条件如下：
初始条件：
（2）
边界条件：
（3）
一般假定函数 ， ， 满足连接条件，即 ， 。
enddo
!-----------!
!-------------------每隔15s输出沿程变化------------------------!
do n=1,nt
t=(n-1)*dt
if (mod(int(t*1000),15000).eq.0) then
write(1,*) 't=',t,'s'
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')
图3 x=10m处u随时间变化图
5.2空间步长及时间步长
通过对差分格式的稳定性分析知，显格式空间步长与时间步长间应满足一定的关系，即式（21）。本文选用的步长如表1所示，分别用显格式和全隐格式进行数值计算，差分网格如图4所示。
表1时间步长及空间步长取值表
时间步长 /s
空间步长 /m
时间步数
计算时长/s
空间步数
0.003
0.1
500000
1500
100
0.3
图4差分网格示意图
5.2求解方法
对于显格式，有 ，若n时间层上的u已知，则可直接求解出（n+1）时间层上的u值。因此，由于给出了初始条件，显格式无需迭代，直接从t=0时刻开始，逐一计算下一个时间层上的u值，便可求解出各时间层上各空间点的u值。
对于全隐格式，由于 与 ， 有联系，不能直接求解，必须联解代数方程组。此时方程组为三对角方程组，可采用追赶法进行求解。
enddo
!-----------!
!----------------------------------求解差分方程-----------------------------!
do n=2,nt
call systri(nx,d,a,c,b,xx)
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)
对一维扩散方程：
（7）
采用显格式差分格式，令 ， ，则
（8）
用Taylor级数展开代入上述差分方程中，则有
（9）
当 ， 时，上式化为
（10）
可见，此差分格式所构成的差分方程与原来的微分方程是相容的，故该显格式为相容格式。