Final Project 有限差分法求解一维对流扩散方程

姚媛媛 2014301020047

摘要

对流扩散方程（convection diffusion equation ）

是一类基本的运动方程，它可以用来对流扩散问题数值计算方法的研究具有重要的理论和实际意义，可用于环境科学、能源开发、流体力学和电子科学等许多领域。我们现在用有限差分的方法来求解一维对流扩散方程并讨论解的精确性。

关键词： 对流扩散方程 理论和实际意义 有限差分 一维

---

背景

问题引入

考虑一维对流扩散方程

$$\left\{ \begin{array}{ll} -au^{\prime\prime}+bu^\prime = 0, & 0 < x < L =1, \\ u(0)=0, u(L)=1. & \end{array} \right.$$
其精确解为
$$u(x)=\frac{1-e^{Rx}}{1-e^R},$$
其中$R=\frac{bL}a$称为$P\acute eclet$数。

构造解法

有限差分格式的构造

令$h=\frac{L}{n+1},x_i=ih,i=0,1,\cdots ,n+1$

（1）则对于中心差分格式：

$$u^\prime=\frac{u_{i+1}-u_{i-1}}{2h},u^{\prime\prime}=\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}$$
且定义$c=\frac{Rh}{2}$,将原方程第一式改写为
$$h^{-2}[-(1-c)u_{i+1}+2u_i-(1+c)u_{i-1}]=0(*)$$
于是对于该方程，我们可以构造线性方程$Ax=b$
其中可得系数矩阵$A=h^{-2}\left[ \begin{array}{ccccc} 2 & -(1-c) & 0 & \cdots & 0 \\ -(1+c) & 2 & -(1-c) & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & \vdots \\ \vdots & \vdots & -(1+c) & 2 & -(1-c) \\ 0 & 0 & \cdots & -(1+c) & 2 \\ \end{array} \right]$
并由上述方程计算出$b=h^{-2}\left[ \begin{array}{r} 0 \\ \vdots \\ 0 \\ 1-c \end{array} \right]$
且$x$表示以元素${u_1,u_2,\cdots,u_n}$构成的列向量，
这样我们就可以通过求解线性方程组得到$u_i$的值，从而拟合曲线逼近精确解

(2)对于迎风格式，分两种情况讨论：

$b>0$时，构造向后差分格式：

$$u^\prime=\frac{u_i-u_{i-1}}{h},u^{\prime\prime}=\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}$$
且定义$c=Rh$,将原方程第一式改写为
$$h^{-2}[-(1+c)u_{i+1}+(2+c)u_i-u_{i-1}]=0(*)$$
于是对于该方程，我们同样可以构造线性方程$Ax=b$
其中可得系数矩阵$A=h^{-2}\left[ \begin{array}{ccccc} 2+c & -1 & 0 & \cdots & 0 \\ -(1+c) & 2+c & -1 & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & \vdots \\ \vdots & \vdots & -(1+c) & 2+c & -1 \\ 0 & 0 & \cdots & -(1+c) & 2+c \\ \end{array} \right]$
并由上述方程计算出$b=h^{-2}\left[ \begin{array}{r} 0 \\ \vdots \\ 0 \\ 1 \end{array} \right]$
且$x$表示以元素${u_1,u_2,\cdots,u_n}$构成的列向量，
这样我们就可以通过求解线性方程组得到$u_i$的值，从而拟合曲线逼近精确解

$b<0$时，构造向前差分格式：

$$u^\prime=\frac{u_{i+1}-u_i}{h},u^{\prime\prime}=\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}$$

且定义$c=Rh$,将原方程第一式改写为

$$h^{-2}[-u_{i+1}+(2-c)u_i-(1-c)u_{i-1}]=0(*)$$

于是对于该方程，我们同样可以构造线性方程$Ax=b$
其中可得系数矩阵$A=h^{-2}\left[ \begin{array}{ccccc} 2-c & -(1-c) & 0 & \cdots & 0 \\ -1 & 2-c & -(1-c) & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & \vdots \\ \vdots & \vdots & -1 & 2-c & -(1-c) \\ 0 & 0 & \cdots & -1 & 2-c \\ \end{array} \right]$
并由方程计算出$b=h^{-2}\left[ \begin{array}{r} 0 \\ \vdots \\ 0 \\ 1-c \end{array} \right]$
且$x$表示以元素${u_1,u_2,\cdots,u_n}$构成的列向量，
这样我们就可以通过求解线性方程组得到$u_i$的值，从而拟合曲线逼近精确解

---

正文

算法描述

$Jacobi$与$Gauss-Seidel$迭代的算法描述
考察非奇异线性代数方程组$Ax=b$,且将$A$拆分成该种形式$A=D-L-U$
其中$A=[a_{ij}],D=diag(a_{11},a_{22},\cdots,a_{nn})$

$$L=\left[ \begin{array}{ccccc} 0 & 0 & 0 & \cdots & 0 \\ -a_{21}& 0 & 0 & \cdots & 0 \\ -a_{31}&-a_{32} & 0 & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & 0 \\ -a_{n,1}& -a_{n,2} & \cdots & -a_{n,n-1}& 0 \\ \end{array} \right]\quad U=\left[ \begin{array}{ccccc} 0 & -a_{12}& -a_{13}& \cdots & -a_{1,n} \\ 0 & 0 & -a_{23}& \cdots & -a_{2,n} \\ 0 &0 & 0 & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & -a_{n-1,n} \\ 0 & 0 & \cdots & 0 & 0 \\ \end{array} \right]$$
此时构造单步线性定常迭代法$x_k=Mx_{k-1}+g$
(1)对于$Jacobi$迭代法，有$M=D^{-1}(L+U),g=D^{-1}b$
若给定初始向量$x_0=\left[ \begin{array}{r} x_1^{0}\\ x_2^{0}\\ \vdots\\ x_n^{0} \end{array} \right]$\
代入方程式的右边，即可算出一个新的向量$x_1$,依次类推，可得精确解的近似解

(2)对于$G-S$迭代法，有$M=(D-L)^{-1}U,g=(D-L)^{-1}b$
注意到$Jacobi$迭代中各分量的计算顺序是没有关系的，
但是$G-S$迭代法是在计算$x_k$的第一个分量用$x_{k-1}$的各个分量计算，但在计算第二个分量时
$x_1^{k}$已经算出，用它代替$x_1^{k-1}$,其他分量仍用$x_i^{k-1}$依次类推
$G-S$迭代法的一个明显好处是在编程时存储量减小
计算该种迭代法，计算分量的次序不可以改变

固定a=1, L=1, 选取不同的b，分别用两种差分格式求解对流扩散方程，然后任选一种迭代格式求解线性方程组，以绘图的方式给出数值结果，并解释相应数值现象。
其精确解为：

结果

(1)取$b=-50$,此时将中心差分格式与向前差分格式对比
可知中心差分格式的拟合程度优于向前差分格式
中心差分格式和向前差分格式迭代次数分别为1143次和3691次,说明中心差分格式收敛更快

(2)取$b=-25$,此时将中心差分格式与向前差分格式对比
可知中心差分格式的拟合程度优于向前差分格式
中心差分格式和向前差分格式迭代次数分别为3111次和9029次,说明依然是中心差分格式收敛更快

(3)取$b=25$,此时将中心差分格式与向后差分格式对比
可知中心差分格式的拟合程度优于向后差分格式
中心差分格式和向前差分格式迭代次数分别为1622次和6308次,说明中心差分格式收敛更快

(4)取$b=50$,此时将中心差分格式与向后差分格式对比
可知中心差分格式的拟合程度优于向后差分格式
中心差分格式和向前差分格式迭代次数分别为435次和2022次,说明中心差分格式收敛更快

(5)取$b=0$,此时将中心差分格式与精确解对比
可知中心差分格式的拟合程度较好，此时中心差分格式可写为$2*u_i=u_{i-1}+u_{i+1}$，即等差数列
中心差分格式迭代次数分别为32914次

(6)我们可以看出中心差分格式都拟合地较好，然而事实上，当$h>\frac{2}{R},\sigma =\frac{1+c}{1-c}<0,u_i$会在零点附近产生震荡，而精确解在$[0,1]$上恒正且单调上升。此时近似解变得不精确，是格式本身带来了震荡，为了回避这一点，我们应选取足够小的$h$以保证$0<c<1$，即$h<\frac{2}{R}$，且当$b$越大就应$h$越小，以减小误差。
如取$b=-50,n=21$,则可发现震荡图像

(7)此时改变矩阵维数，情况如下：

程序代码：CODE

---

结论

从图中我们可以看出中心差分格式都拟合地较好，其精确度比向后差分格式要高。

---

致谢

感谢蔡浩老师一个学期以来的辛勤教导，感觉接触了很多很牛逼的东西，可是因为本身基础有限，之前没有接触过计算这一块的东西，所以到期末会的东西也很简陋。虽然以后应该是不会专门朝计算物理这个方向发展，但即使不做这些，现在学的也对以后很有用，也是很有必要掌握的。
虽然说期末程序得自己写，但是动笔之前还是看了好多学长学姐的期末作业的，五花八门，各式各样。然后可能也有点启示吧，在此感谢各位大佬们。
最后还要感谢数院的颜值担当，最美丽的sys帮我解决了数学上的问题，要不然这些让人头疼的矩阵根本就看不懂什么意思。