最近做耦合问题又需要用到有限差分,就把这个放在草稿箱里的有限差分整理一遍发出来就当复习了!
1.有限差分法的引入与介绍
如果用一句话来概括有限差分法,那就是利用用网格节点逼近导数并建立代数方程组来求解各个网格上节点的值。在系列一里面我们知道要把一个无限、连续的问题通过区域的划分变成有限、离散的问题。在有限差分里的区域划分是等距的网格。
有限差分:利用用网格节点逼近导数,并建立有限个未知数的代数方程组来求解各个网格上节点的值。
上面这句话很显然就产生了两个问题,如何用网格节点逼近导数?(可以用高数里面学的泰勒展示)逼近哪里的导数?(这里面的区别就产生了向前向后和CN)。解决了上述两个问题,也就产生了不同的差分格式。
1.离散区域并给出差分格式
获得差分格式其实并不是一个很复杂的过程,随便几种差商格式组合、泰勒展开都能得到很多的差分格式,但有限差分最核心的地方并不在于差分格式构造,而是在于分析得到的差分格式靠谱不靠谱。
而这个分析格式是否靠谱,这就像兄弟们找对象,对方能不能逼近自己的理想型,这个叫收敛性,现实的遇到的那个人(数值解)和自己的理想型(方程的真解)越接近当然越好;在一起后遇到小矛盾(对应过来是小扰动或者截断、舍入误差)两个人的处理方式和后续的发展,这个叫稳定性,不管多小的矛盾都要闹成大矛盾,这种就是不稳定的,当然得淘汰掉ಠ_ಠ。
收敛性这个没有什么好说的,在本系列第一篇文章中有说到。至于稳定性,电脑计算小数时四舍五入造成的误差(舍入误差)一定存在,就要考虑他在逐层求解的时候是否会变得越来越大。常用的分析方法有矩阵法和fourier方法。
相容性是在得到一个点的值后,他和初始定义的值能不能相等,这个是相容性,一般是在第四步代数方程组求解完在验证。
以上一起就是有限差分法最核心的理论分析:
2.差分解的唯一性、收敛性计稳定性的讨论
接下来,这一步就是在前面的离散里面得到了一堆离散的关系表达式,如果这个关系表达式是显式的,而第一层(初边值条件)上的每一个离散点是已知的,这样可以用第一层离散点值来直接表示出第二层上所有的离散点的值。依次这样算下去,所有层上的离散点都能被求出来。当网格足够密的时候,也就是离散点足够多的时候,可以将其近似看作求出了 U U U在区域内的解。
对隐式差分格式的计算,一般来说是一个三对角矩阵线性代数方程组的求解。这个问题就要从差分理论的战场转移到线性代数求解算法和计算复杂度优化上了:在数学上可以用追赶法和SOR方法来求解优化,在计算机上可以用稀疏矩阵的格式来存储来优化(因为三对角阵的0元素比较多)
这就是理论分析过后进行实际求解并且和计算机结合起来的一步:
3.求解代数方程组
最后一步当然就是通过程序的结果看看他是否和我们当时理论所计算的误差精度、稳定性是否一样啦
4.验证程序结果和理论结果的一致性
2.几个常见的差分格式
2.1 一维问题
1.对求解域做网格剖分
我们将区间[a,b]划分为等距的n个部分,对应的有n+1个节点,取这些节点为:
x 0 , x 1 . . . x i . . . x n x_0,x_1...x_i...x_n x0,x1...xi...xn
相邻节点间距为 h = b − a n h=\frac{b-a}{n} h=nb−a,记为步长。
2.几种常见的差分格式
差分格式:用离散网格点的组合(这里也就是差分或者差商)代替方程中的偏导数的格式
在一维情况,取等距模板点 x i x_i xi、 x i − 1 x_{i-1} xi−1、 x i + 1 x_{i+1} xi+1,两点间距都为 h h h,常见的差分格式有:
- 向前差分格式:即两个离散点的差商是用来代替前面那个点的导数
u x = x i + 1 − x i h u_{x}=\frac{x_{i+1}-x_{i}}{h} ux=hxi+1−xi - 向后差分格式:即两个离散点的差商是用来代替后面那个点的导数
u x = x i − x i − 1 h u_{x}=\frac{x_{i}-x_{i-1}}{h} ux=hxi−xi−1 - 中心差分格式 :这时候要三个离散点,最边上两个点的差商来代替中间那个点的导数。
u x = x i + 1 − x i − 1 2 h u_{x}=\frac{x_{i+1}-x_{i-1}}{2h} ux=2hxi+1−xi−1
这时候也能刻画二阶导数,这就是二阶中心差分格式
u x x = x i + 1 + x i − 1 − 2 x i h 2 u_{xx}=\frac{x_{i+1}+x_{i-1}-2x_{i}}{h^2} uxx=h2xi+1+xi−1−2xi
热传导方程中C-N格式(克兰克-尼科尔森方法)必须得拿来介绍一下。因为是在热传导方程中,所研究的对象 u u u 就是描述一个区域内的温度如何随时间变化,所以他含有对时间的偏导数。
CN格式的核心思想简单理解就是在时间层 t = n + 1 / 2 t = n+1/2 t=n+1/2处的导数用中心差分格式,在 x j xj xj用二阶中心差分格式。他的思路来源是:既然(j,n+1)(j,n)两点的差商可以代替t = n和t = n+1处的导数,那么能不能用这两个点的差商来代替中间点 t = n + 1 / 2 t = n+1/2 t=n+1/2的导数呢?当然是可以的,并且代替后,因为 t = n + 1 / 2 t = n+1/2 t=n+1/2这个时刻虽然在现实意义本身是存在的,但是他并没有出现在我们离散的网格点里面,我们就不能把它当作存在来处理,所以就用中心差分格式把这个时刻的导数用(j,n+1)(j,n)两点的差商来表示,这就是CN格式的由来。
同理还可以代替其他点比如 t = n + 1 / 12 t = n+1/12 t=n+1/12处的导数,代替不同点的导数截断误差也会有差别。
注1:在二维情况下,根据不同的目标问题,还可以取更多的模板点来构造差分格式以获得更好的数值解,例如椭圆型问题中取六个点的六点对称格式、激波管问题中引进人工粘性项的Lax-Wendroff格式等。
(待更新)
综上,所以可以看出来有限差分的步骤并不难,理论的基础相比有限元法要简单很多,以后的文章再利用具体的代表算例来求解。