(1)在数学中,有限差分方法(Finite Differenc Methods,FDM),是一种微分方程的数值方法,是通过有限差分来近似导数或者偏导数,从而求得微分方法的近似解。在实际中,许多数学问题都很难得到其解析解,所谓解析解就是可以通过给出解的具体函数形式,从解的表达式中就可以算出任何自变量对应值;而数值解是在特定条件下通过近似计算得出来的一个数值。 举个例子: 一般一元二次方程 a x^2 + b x + c == 0 (a ≠ 0) 的解析解为 x = (- b ±√(b^2 - 4 a c)) / (2 a) 数值解只能由具体的系数确定。 (2)导数的有限差分近似由泰勒展开的推算 假设f(x)在x0处各阶可导,则根据定理,其在x0+h的泰勒展开式为: {\displaystyle f(x_{0}+h)=f(x_{0})+{\frac {f'(x_{0})}{1!}}h+{\frac {f^{(2)}(x_{0})}{2!}}h^{2}+\cdots +{\frac {f^{(n)}(x_{0})}{n!}}h^{n}+R_{n}(x),} 其中 n !表示是 n 的 階乘 , R n ( x )為餘數,表示泰勒多項式和原函數之間的差。可以推導函數 f 一階導數的近似值: 设定x 0 =a,可得: {\displaystyle f(a+h)=f(a)+f'(a)h+R_{1}(x),} 除以 h 可得: {\displaystyle {f(a+h) \over h}={f(a) \over h}+f'(a)+{R_{1}(x) \over h}} 求解f'(a): {\displaystyle f'(a)={f(a+h)-f(a) \over h}-{R_{1}(x) \over h}} R1(x)是x(x-a时)的高阶无穷小,假设R1(x)很小并可以忽略, 可以将f的一阶导数近似为: (一阶近似) 定理 : 设 n 是一个 正整数 。如果定义在一个包含 a 的 区间 上的函数 f 在 a 点处 n +1 次 可导 ,那么对于这个区间上的任意 x ,都有: {\displaystyle f(x)=f(a)+{\frac {f'(a)}{1!}}(x-a)+{\frac {f^{(2)}(a)}{2!}}(x-a)^{2}+\cdots +{\frac {f^{(n)}(a)}{n!}}(x-a)^{n}+R_{n}(x).} 其中的多项式称为函数在 a 处的 泰勒展开式 ,剩余的 {\displaystyle R_{n}(x)} 是泰勒公式的余项,是 {\displaystyle (x-a)^{n}} 的高阶无穷小。 由于在近似过程中,我们截断了Rn(x),因此数值差方法本来存在误差,这种由于截断所引起的误差就是结尾误差。 差分的方法有不同分类,如前向,后向和中心差分: 前面两种是对偏微分方程一阶近似,中心差分式二阶近似。 差分还可以分为,explicit and implicit difference 显式方法 热传导方程最常用显式方法的 模版 利用在时间 {\displaystyle t_{n}} 的前向差分,以及在位置 {\displaystyle x_{j}} 的二阶中央差分( FTCS 格式 ),可以得到以下的迭代方程: {\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}.\,} 这是用求解一维导热传导方程的 显式方法 。 可以用以下的式子求解 {\displaystyle u_{j}^{n+1}} {\displaystyle u_{j}^{n+1}=(1-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}} 其中 {\displaystyle r=k/h^{2}.} 因此配合此迭代关系式,已知在时间 n 的数值,可以求得在时间 n +1的数值。 {\displaystyle u_{0}^{n}} 及 {\displaystyle u_{J}^{n}} 的数值可以用边界条件代入,在此例中为0。 此显式方法在 {\displaystyle r\leq 1/2} 时,为 数值稳定 且 收敛 。其数值误差和时间间隔成正比,和位置间隔的平方成正比: {\displaystyle \Delta u=O(k)+O(h^{2})\,} 隐式方法 隐式方法的模版 若使用时间 {\displaystyle t_{n+1}} 的后向差分,及位置 {\displaystyle x_{j}} 的二阶中央差分(BTCS 格式),可以得到以下的迭代方程: {\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}}.\,} 这是用求解一维导热传导方程的 隐式方法 。 在求解线性联立方程后可以得到 {\displaystyle u_{j}^{n+1}} : {\displaystyle (1+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=u_{j}^{n}} 此方法不论 {\displaystyle r} 的大小,都数值稳定且收敛,但在计算量会较显式方法要大,因为每前进一个时间间隔,就需要求解一个联立的数值方程组。其数值误差和时间间隔成正比,和位置间隔的平方成正比: {\displaystyle \Delta u=O(k)+O(h^{2})\,} (3) 舍入误差是由于计算机数字限制(四舍五入)引起的 由于计算机的 字长 有限,进行数值计算的过程中,对计算得到的中间结果数据要使用“四舍五入”或其他规则取 近似值 ,因而使计算过程有 误差 。这种误差称为舍入误差。 在python中输入0.1+0.2,得到的结果是0.30000000000000004,这是舍入误差带来的结果。 There are two approximations in numerical modelling. One needs to ask the questions: How well is the natural system modelled by the differential equations?, and, How well is the solution to the differential equations represented by the computational algorithm?. I n the analysis here more attention is paid to the second question. The first question can only be answered by studying the behaviour of the natural system and comparing i t to the equations applied to it. Therefore i t will be assumed here that the differential equations approximate the system we1 I despite the fact that i t has been noticed that this is not necessarily the case. Abbott (1974) noticed that a difference scheme considerably different from the differential equations used to describe a system, can yield more accurate results than a difference scheme similar to the differential equations when compared with experimental results. There are three possible sources of error associated with finite difference solutions to partial differential equations. I t is important that one understands these sources, their consequence if not control led, and means for controlling them. These three sources of error are: truncation, discretization, and round-off. Truncation error occurs when a derivative is replaced with a finite difference quotient; discretization error is due to the replacement of a continuous model (function) with a discrete model; and round-off error is essentially machine error in that the algebraic finite difference equations are not always solved exact l y . For finite difference solutions to be accurate, they must be con- sistent and stable. Consistency simply means that the truncation errors tend to zero as Ax and At - 0, i.e., as Ax and At - 0 the finite difference equation becomes the original differential equation. This is examined i n the follcwing paragraphs. Stability implies the controlled growth of round-off error. Stability considerations apply principally to explicit schemes to be discussed later. Any numerical scheme that allows the growth of error, eventually swamping the true solution, is unstable. References (1) https://zh.wikipedia.org/zh-cn/%E6%9C%89%E9%99%90%E5%B7%AE%E5%88%86%E6%B3%95 (2) https://zh.wikipedia.org/wiki/%E6%B3%B0%E5%8B%92%E5%85%AC%E5%BC%8F (3) Stephenson, D., kinematic hydrology and modelling. 1986. (4) Chow, E.A., Applied hydrology. 1988.
数值差分方法一般分为显式差分与隐式差分,对二者的区别我查阅了相关的资料。我后面发现关于二者区分不是有一些差异。如Wiki: Explicit methods calculate the state of a system at a later time from the state of the system at the current time, while implicit methods find a solution by solving an equation involving both the current state of the system and the later one. Mathematically, if {\displaystyle Y(t)} is the current system state and {\displaystyle Y(t+\Delta t)} is the state at the later time ( {\displaystyle \Delta t} is a small time step), then, for an explicit method {\displaystyle Y(t+\Delta t)=F(Y(t))\,} while for an implicit method one solves an equation {\displaystyle G{\Big (}Y(t),Y(t+\Delta t){\Big )}=0\qquad (1)\,} to find {\displaystyle Y(t+\Delta t).} 下面几个链接是其它的一些讨论: http://earthscience.stackexchange.com/questions/2253/what-are-the-differences-between-implicit-and-explicit-methods-when-applied-to-m http://www.cfd-online.com/Forums/main/5982-implicit-vs-explicit-method.html 这几个讨论,认为如果下一时刻状态值,能通过当前时刻值进行直接计算,则表示显式差分,如果下一时刻状态值,不能通过当前时刻值进行直接计算,而必须求解(1)式来计算,则表示隐式差分。 这个定义似乎很简单,也容易理解,但现在有一个情况是,式(1)如果是线性的,并且通过式(1),我们能反解出Y(t+dt),使其能够写出显示差分的格式。举一个例子吧, du/dt = f(u) where u is a vector and f is a known function. Then we may approximate the time derivative over a time step h as (implicit) u(n+1) = u(n) + h.f(u(n+1)) if f(u) = Au where A is a constant, then u(n+1)=u(n)+Ahu(n+1),then, u(n+1)=u(n)/(1-Ah) 这种情况下,差分方法是显式还是隐式呢? 本质上我们也求解方程(1), 但是我们得到了显示格式。 后面我又查了相关的资料,如: https://www.flow3d.com/home/resources/cfd-101/numerical-issues/implicit-versus-explicit-numerical-methods 这里认为: When a direct computation of the dependent variables can be made in terms of known quantities, the computation is said to be explicit. When the dependent variables are defined by coupled sets of equations, and either a matrix or iterative technique is needed to obtain the solution, the numerical method is said to be implicit. 根据这个理解,如果待求的值,能用已知的(已经计算的)值直接计算出来,则表示 显式差分,相反,则表示隐式差分。 如下面的例子,我们计算某一位置某一时刻的值时,不仅用到了前一时刻的信息,也用到了当前时刻的信息,最后通反求解的方式,得到计算值,然而此类差分仍然属于显式差分。 另外一个定义: The idea behind the explicit method is to be able to obtain values such as y i+1 = f(xi, y i ), y i+2 = f(x i ,x i+1 ,y i ,y i+1 ) , etc. In other words, your solution proceeds by solving explicitly for a new unknown value in the solution array, given all previous values in the array. On the other hand, implicit methods imply the simultaneous solution of n linear algebraic equations that provide, at once, the elements of the solution array. 对于隐式差分,我有个例子 在上面的式子中,我们并不能直接求解出yi,这时候需要更多的方程来求解,这些方程可以写成矩阵来求解。 关于两者的Advantage与Disadvantage,应该比较清楚,前者速率快,但是本质上是conditionally stable,当dt较大,精度就会受到限制;而后者是unconditionally stable,但是计算效率低。 1-s2.0-S0096300304004680-main.pdf role_959843.pdf tcs-30_web.pdf 今天查看另外一个链接: http://echo2.epfl.ch/VICAIRE/ ,相关信息如下: In finite-difference methods, calculations are performed on a grid placed over the interesting flow domain into the x-t plane. This computational grid is defined by some equal or variable space and time steps, Δ x and Δ t , respectively. A network of discrete points is thus obtained and the flow variables Q, z (or Q, h; V, z; V, h ) are derived only for this finite number of grid points. Figure 6.2 shows a typical computational grid. The spatial positions of the grid points are denoted by index j and the time moments by index i . Figure 6.2. Finite-difference computational grid. As principle, finite-difference methods transform governing partial differential equations into a set of algebraic (linear or nonlinear) finite-difference equations, which are solved to allow the values of flow variables in a grid point or in all grid points on a time line. These finite-difference equations are derived by approximating the time and space derivatives with some finite-difference expressions. Not only the derivatives, but also the other terms into the Saint-Vénant equations must be defined in a certain manner. Concerning the space derivative of a continuous function f(x,t) at time moment t i and space position x j on the grid in Figure 6.2, this can be approximated as: a forward difference approximation; a backward difference relationship; a central difference expression, (6.29) in which f j i represents the value of f at grid point ( x j ,t i ). In an analogous manner, the time derivative may be defined in several different ways as, for example: ; (6.30) Usually, the value of the function f(x,t) is accepted at grid point ( x j ,t i ) as f j i , but some different approximations can also be decided. A finite difference method must employ a certain type of finite-difference scheme. These schemes are grouped in two major classes: explicit and implicit finite-difference schemes. Explicit schemes are those in which the flow variables at any point j and time level i +1 may be computed using only known data at a few adjacent points on the time line i . These schemes do not lead to a system of algebraic equations, but rather to only two finite difference equations for each grid point ( x j ,t i+1 ). By solving the two equations the unknown values of flow variables are obtained and then, the computation proceeds to the next grid point along the time line i +1. In implicit schemes, finite-difference expressions used to approximate the space and/or time derivatives at grid point ( x j , t i+1 ) include the unknown values of flow variables at a few adjacent points on the time line i +1. Consequently, a system of algebraic equations is produced for a given time line i +1 and by solving this system, all unknown values are simultaneously determined at time level i +1. Replacing the continuous original problem with an integration over a discrete computational grid introduces numerical errors into the results. A finite-difference scheme is stable if such errors are not amplified during computation from one time level to the next. The numerical stability depends on the size of the time and space steps and on some flow characteristics. The Courant condition: (6.31) is a necessary but insufficient condition for stability of an explicit scheme. Therefore, any explicit scheme is conditionally stable, the Courant condition requiring to work with small time steps as compared with the physical phenomena evolution. Despite their computational simplicity, the explicit methods are seldom used in river modelling for reasons of this stability restriction. On the other hand, the implicit schemes may appear more complicated, but can generally be made unconditionally stable for large computation steps and with little loss of accuracy. Another distinction among schemes belonging to the same class is related to the way in which the non derivative terms (such as S f (Q,h) , B(h) , A(h) etc.) are discretized. Because these terms are functions of dependent variables, their treatment induces a linear or nonlinear feature of the finite-difference equations. A lot of finite-difference schemes will be presented in more details within the next sections. Numerical solution of kinematic wave model First form (6.21) of the kinematic wave model is used to illustrate an explicit linear scheme on a finite-difference cell as show in Figure 6.4. A backward finite-difference is used to approximate the space derivative, while the time derivative is usually expressed at the same x j+1 spatial position. In order to obtain a linear scheme, the nonderivative term α β Q β -1 must be evaluated for a known Q value, which here is accepted as mean between the two diagonal values Q j i+1 and Q i j+1 . Figure 6.4. Finite-difference cell for linear kinematic wave solution. Consequently, the finite-difference form of equation 6.21 is: (6.40) and the unknown Q j+1 i+1 results as: (6.41) Starting with the inflow hydrograph value Q at the time level t i+1 , the computation sequentially proceeds from upstream towards downstream grid points, along the current time line i +1. However, this scheme supposes that within the coefficient α , the wetted perimeter P remains constant. On the other hand, by using Q instead of A into the time derivative ∂ A/ ∂ t , relative computational errors are decreased. Indeed, taking the logarithm of equation 6.20, i.e. , and differentiating, one obtains: , where β is 0.6 for the Manning's equation. It follows that the estimation error in Q would be magnified by 1/0.6 ≈ 1.67 if the cross-sectional area A is used as dependent variable. 来自另外一本书:Chow, E.A., Applied hydrology. 1988. A finite-difference method may employ either an explicit scheme or an implicit scheme for solution. The main difference between the two is that in the explicit method, the unknown values are solved sequentially along a time line from one distance point to the next, while in the implicit method the unknown valueson a given time line are all determined simultaneously. The explicit method is simpler but can be unstable, which means that small values of Ax and A^ are required for convergence of the numerical procedure. 综合上面的分析,对于PDE,如果在某个待求得时间点上,各位置(from upper to down)的值能够被依次直接求解的情况,认为是显示差分;而如果 某个待求得时间点上,各位置的值必须同时联立起来求解,应该是隐式差分。