Direct versus iterative methods for forward-backward diffusion equations. Numerical comparisons

Direct versus iterative methods for forward-backward diffusion equations. Numerical comparisons on a particular transport kinetic model. ´ Oscar L´opez Pouso∗ Nizomjon Jumaniyazov† October 5, 2020 Abstract By far, the standard implementation of finite difference schemes for forward-backward partial differential equations consists in employing an iterative method. This paper collects a series of numerical results which demonstrate that a direct implementation can reduce the computing time. An effective way of choosing the seed for the iterative method naturally arises. MSC 2020: Primary 65M06, 65Z05; Secondary 35K65, 35Q84, 78M20, 78A35. Keywords: forward-backward partial differential equation, two-way partial differential equation, finite difference method, Fokker-Planck equation, direct method, iterative method. 1 Introduction This work has been motivated by reference [11], which proposes a couple of finite difference schemes for solving certain problem of interest to the nuclear engineering community. To specify, the problem is the one defined on Q= [−1,1] ×[Zini, Zfin] by the Fokker-Planck equation µ∂ψ ∂z +αψ −σ∂ ∂µ [D(µ)∂ψ ∂µ ]=Wfor (µ, z)∈Q, (1) being D(µ) = 1 −µ2, and the incoming flux conditions ψ(µ, Zini) = f(µ) for µ∈(0,1],(2) ψ(µ, Zfin) = g(µ) for µ∈[−1,0),(3) ∗Dept. of Applied Mathematics, Faculty of Mathematics, University of Santiago de Compostela. C/ Lope G´omez de Marzoa s/n, Campus Vida, 15782 Santiago de Compostela (A Coru˜na, Spain). Email: [email protected]. †Dept. of Natural Sciences, Urgench branch of Tashkent University of Information Technologies named after Muhammad al-Khwarizmi (Uzbekistan). Email: nizomjon [email protected]. 1 where α,σand Ware known functions of (µ, z) with α≥0 and σ > 0, and fand gare known functions of µ. The following paragraph, excerpted from [11], is intended for those readers interested in theoretical considerations: “The fact that the µ= 0 extreme is open in conditions (2) and (3) as well as the absence of boundary conditions at |µ|= 1 can be explained by physical considerations [...] From the mathematical viewpoint, there are also analytical reasons that support the well-posedness of Problem (1)–(3) in many particular cases; the interested reader might like to consult references [2], [3], [5], [6], [10], [12], and [13]. In particular, the absence of boundary conditions at |µ|= 1 is related to the degeneracy of the inner diffusion coefficient 1 −µ2at |µ|= 1.” Again in [11] it is explained that the schemes defined therein admit of two different implementations, one of them which is of iterative type and another one which is not, and which we shall refer to as direct. The precise meaning that both concepts, iterative and direct, have in the present context will be clarified in subsequent paragraphs. Then, the direct implementation is chosen, and no comparison with the other approach in terms of computing time is made. From a general perspective, this equation belongs to the category of the forward-backward or two-way partial differential equations (PDEs). It is remarkable the fact that, before [11], the existing literature on the numerical resolution of this kind of problems by means of finite difference methods has given prominence to propose and analyze iterative algorithms (see [7], [8], [14], [15], [16]). For this reason it is natural to investigate whether the iterative approach results in a faster method or not; the major contribution of this paper is to provide a body of numerical evidence that in fact it is slower in most cases. With the aim of illustrating the matter in question, let us consider the simple forward-backward heat equation x∂u ∂t −∂2u ∂x2=f(x, t),(x, t)∈[−1,1] ×[Tini, Tfin].(4) The PDE (4) is parabolic when x > 0 and backward parabolic when x < 0. Consequently, leaving apart the necessity of imposing some type of boundary conditions at |x|= 1, it demands to split the spatial domain [−1,1] into two parts, (0,1] and [−1,0), in order to impose an initial condition on (0,1] and a final condition on [−1,0) (see [16]). From a numerical “finite differences” standpoint, it is clear that this problem cannot be solved with a marching method by means of one single sweep, as it is the case for the heat equation ∂u ∂t −∂2u ∂x2=f(x, t); indeed, for the equation (4) it is not possible to advance starting at time t=Tini because the initial condition only holds on (0,1], and it is not possible either to go back starting at time t=Tfin because the final condition only holds on [−1,0). As said above, the standard treatment is of iterative type, and consists in: 1. Providing a seed, that is to say, a guess of the solution on the interface {0} × [Tini, Tfin], 2. Solving on [0,1] ×[Tini, Tfin], advancing in time, 3. Solving on [−1,0] ×[Tini, Tfin], going back in time, 4. Updating the value on the interface, and 2 5. Iterating steps 2, 3 and 4 until convergence is achieved. It is not being claimed that all kind of forward-backward PDEs can be solved in this way. In fact, there are more complex examples that must viewed from some other perspective (see [1]), but the problem (1)–(3) admits of a similar treatment. As an alternative to the method above, one has what it is called in this paper the direct approach, consisting in solving the problem simultaneously on the whole time-space domain, treating time tas another spatial variable, as if one were solving, in the elliptic manner, the Dirichlet problem for the Poisson equation in two dimensions. The principle of the iterative approach is that it is much faster to march in time than to solve on the whole time-space domain by means of the direct approach. In other words, it is much faster to solve Nlinear square systems of order Ithan one single linear square system of order I×N. However, since in this case marching in time implies performing iterations, the following question arises: is it possible that the number of iterations caused the iterative method to be slower than the direct method? This is a point of interest that, to our knowledge, has not been explored yet. As said above, this paper focuses on showing that the direct approach is generally advantageous over the iterative approach when solving the problem (1)–(3) by means of finite differences. It is important to say that the scope of this assertion goes beyond the example studied, in the sense that the same statement is true for other forward-backward diffusion problems, such as the boundary value problem studied in [7], the forward-backward heat equation (4) or, more generally, as it stands in [16] when a=a(x). To the authors knowledge there is no reference, apart from [11], where the direct method, as it is being understood here, be used. The paper is organized as follows: Section 2 establishes some basic hypotheses and notations. Sections 3 and 4 include, respectively, descriptions of the direct and iterative algorithms. Section 5 collects the numerical results and their analysis. Some further explanations are given in Section 6, and finally Section 7 contains the conclusions. 2 Prefatory comments It is clear that in Problem (1)–(3) zis the time-like variable and µis the space-like one. Then, Equations (2) and (3) can be interpreted, respectively, as an initial and a final condition. In respect of hypotheses on the data functions, we shall assume that α, σ and Ware continuous on Q, α ≥0, σ > 0,(5) fis continuous on [0,1], g is continuous on [−1,0].(6) Recall that Q= [−1,1] ×[Zini, Zfin]. Let {(µi, zn) : i∈ {1,...,I}, n ∈ {1, . . . , N}} be a mesh of Qobtained as the Cartesian product of two uniform meshes of [−1,1] and [Zini, Zfin], and think of using a finite difference scheme over this mesh. Let h=2 I−1 be the distance between µ-nodes and let k=Zfin −Zini N−1be the distance between z-nodes. Among all possible finite difference schemes, we choose the odd scheme described in [11], which is the best option out of the two schemes studied 3 therein. Consequently, I(the number of µ-nodes) is assumed to be odd, and µi⋆= 0 if i⋆=I+1 2. The following notations will be also employed: Di=D(µi), Di±1 2= D(µi±h 2); αn i=α(µi, zn), σn i=σ(µi, zn), Wn i=W(µi, zn); fi=f(µi), gi=g(µi); ψn i≈ψ(µi, zn), understanding that ψis the exact solution of Problem (1)–(3). 3 Description of the direct algorithm As said in Section 2, the finite difference scheme we choose is the so-called odd scheme of [11], which exhibits order 2 with respect to each of the two variables when the exact solution is regular on the compact set Q. It reads as follows (from Equation (7) to Equation (12)): •For (i, n)∈ {1} × {1, . . . , N −1}, (−µ1 k+αn 1 2+σn 1D2 2h2)ψn 1+(−σn 1D3 8h2)ψn 2+ +(−σn 1D2 2h2)ψn 3+(σn 1D3 8h2)ψn 4+ +(µ1 k+αn+1 1 2+σn+1 1D2 2h2)ψn+1 1+ +(−σn+1 1D3 8h2)ψn+1 2+(−σn+1 1D2 2h2)ψn+1 3+ +(σn+1 1D3 8h2)ψn+1 4=Wn 1+Wn+1 1 2.(7) •For (i, n)∈({2, . . . , i⋆−1} ∪ {i⋆+ 1, . . . , I −1})× {1,...,N −1}, (−σn iDi−1 2 2h2)ψn i−1+ + −µi k+αn i 2+ σn i(Di−1 2+Di+1 2) 2h2 ψn i+ +(−σn iDi+1 2 2h2)ψn i+1 +(−σn+1 iDi−1 2 2h2)ψn+1 i−1+ +  µi k+αn+1 i 2+ σn+1 i(Di−1 2+Di+1 2) 2h2 ψn+1 i+ +(−σn+1 iDi+1 2 2h2)ψn+1 i+1 =Wn i+Wn+1 i 2.(8) •For (i, n)∈ {i⋆} × {2, . . . , N −1}, (−σn i⋆ h2)ψn i⋆−1+(αn i⋆+2σn i⋆ h2)ψn i⋆+(−σn i⋆ h2)ψn i⋆+1 = =Wn i⋆.(9) 4 •For (i, n)∈ {I} × {1, . . . , N −1}, (σn IDI−2 8h2)ψn I−3+(−σn IDI−1 2h2)ψn I−2+ +(−σn IDI−2 8h2)ψn I−1+(−µI k+αn I 2+σn IDI−1 2h2)ψn I+ +(σn+1 IDI−2 8h2)ψn+1 I−3+(−σn+1 IDI−1 2h2)ψn+1 I−2+ +(−σn+1 IDI−2 8h2)ψn+1 I−1+ +(µI k+αn+1 I 2+σn+1 IDI−1 2h2)ψn+1 I=Wn I+Wn+1 I 2.(10) •For (i, n)∈ {i⋆,...,I} × {1}, ψ1 i=fi.(11) •For (i, n)∈ {1,...,i⋆} × {N}, ψN i=gi.(12) The direct algorithm, which is well defined whenever I≥5, Iodd, and N≥3, consists in solving the previous linear system of I×Nequations for the I×Nunkowns ψn i. We would like to mention that, at the time the paper [11] was written, we were not aware of the beautiful work [7] by Joel N. Franklin and Eugene R. Rodemich, where the main ideas leading to the odd scheme were already stated (specifically, within the subsection 6.2 on the balanced difference method). Let these lines serve as a tribute to them. It must be noticed that the idea of “direct method” which concerns us is different from the one used in [7]. 4 Description of the iterative algorithm Consider the following subsets of Q:Q−= [−1,0) ×[Zini, Zfin], Q+= (0,1] ×[Zini, Zfin], and Q0={0} × [Zini, Zfin] (see Figure 1). [ ] f(µ) z=Zfin z=Zini µ=−1µ= 1 µ= 0 Q+ g(µ) Q− ) ( Figure 1: Domain Q=Q−∪Q0∪Q+, with indication of the incoming flux boundary conditions. The dotted vertical line represents the subset Q0. 5 STEP 1. Computation of values on Q+at iteration q+ 1: For a fixed number q∈N∪ {0}, assume that (ψn i⋆)[q]for n∈ {2,...,N −1}(15) are known, and obtain all values on Q+at iteration q+ 1, i. e., (ψn i)[q+1] for (i, n)∈ {i⋆+ 1, . . . , I} × {1, . . . , N},(16) by solving Equations (8, for i∈ {i⋆+ 1, . . . , I −1}), (10), and (11), with the aid of the boundary values (15). This is done by means of a one-step forward marching procedure starting from the known values at z=Zini given by Equation (11). STEP 2. Computation of values on Q−at iteration q+ 1: Obtain all values on Q−at the same iteration q+ 1, i. e., (ψn i)[q+1] for (i, n)∈ {1, . . . , i⋆−1} × {1, . . . , N},(17) by solving Equations (7), (8, for i∈ {2, . . . , i⋆−1}), and (12), with the aid of the same boundary values (15). This is done by means of a one-step backward marching procedure starting from the known values at z=Zfin given by Equation (12). STEP 3. Update (computation of values on Q0at iteration q+ 1): Update of the values on Q0can be done by using Equation (9). That is to say, the updated values (ψn i⋆)[q+1] for n∈ {2,...,N −1}(18) can follow from (−σn i⋆ h2)(ψn i⋆−1)[q+1] +(αn i⋆+2σn i⋆ h2)(ψn i⋆)[q+1] + +(−σn i⋆ h2)(ψn i⋆+1)[q+1] =Wn i⋆.(19) However, it is possible to accelerate convergence by over-relaxing the update, in the way we explain now: firstly, one computes (˜ ψn i⋆)[q+1] for n∈ {2,...,N −1}(20) from (−σn i⋆ h2)(ψn i⋆−1)[q+1] +(αn i⋆+2σn i⋆ h2)(˜ ψn i⋆)[q+1] + +(−σn i⋆ h2)(ψn i⋆+1)[q+1] =Wn i⋆,(21) and, finally, the updated values are computed from (ψn i⋆)[q+1] =ω(˜ ψn i⋆)[q+1] + (1 −ω) (ψn i⋆)[q](22) for n∈ {2, . . . , N −1}, where ω∈Ris a given relaxation parameter. Notice that ω= 1 means that no relaxation is being deployed. The value of ωcannot be 0, otherwise no update is taking place, but being different from 0 is not guarantee of convergence. Among those values 7 of ωoffering convergence, the optimal one should be employed, but as of today the issue of finding a closed expression for this optimum (or for an estimation of it) has not been investigated; accordingly, the values of ωused in Section 5 when reporting the numerical results have been found by performing several trials, looking for values that reduce the number of iterations in a significant way with respect to the case ω= 1. The idea of relaxing the update in forward-backward diffusion problems has already been used in [15] and [16]. STEP 4. Checking convergence: Let us define the vectors u[q],u[q+1] ∈RN−2as follows: u[q] j=(ψj+1 i⋆)[q], u[q+1] j=(ψj+1 i⋆)[q+1] (23) for j∈ {1, . . . , N −2}. Given a small real number ε > 0, we say that the algorithm has converged with ε-tolerance if ∥u[q+1] −u[q]∥ 1 + ∥u[q+1]∥≤ε, (24) where ∥ · ∥ stands for the Euclidean norm in RN−2. The quotient in Equation (24) will be referred to as the residual. In case that condition (24) is not satisfied, one sets q=q+ 1 and goes back to STEP 1. If on the contrary condition (24) holds, then the process finishes and the last computed values are taken as the ultimate numerical solution: ψn i= (ψn i)[q+1] for (i, n)∈ {1, . . . , I} × {1, . . . , N}.(25) Notice that both the solution (25) and the number of iterations needed for getting it depend on the seed, on the relaxation parameter ω, on the tolerance parameter ε, and on the norm used in the stopping criterium (24). In what follows, an iteration will be the process which consists in performing steps 1, 2, 3, and 4. 5 Numerical results In this section, Zini = 0, Zfin = 1, and Eabs(Q) = maxQ|ψgrid −ψ|, where ψgrid is representing the approximate solution and the maximum is taken over the set of all nodes. Computations have been performed by running Matlab R ,R2015a, on a personal computer with an Intel R  CoreTM i7-4790 @ 3.60GHz processor. Matlab R automatically chooses, via the single character “\”, direct methods for solving the linear systems involved: LAPACK for the linear systems of the iterative method, and UMFPACK for the only linear system of the direct method. More information about these routines can be found in the reference [4]. All matrices involved in the resolution process are sparse and banded. In the case of the direct method, our code vectorizes the definition of the matrix in order to avoid loops, and, as it is natural, the ordering of the equations and unknowns is chosen so as to produce small bandwidth. 8 We remark that the matrices needed in the iterative algorithm are much easier to deal with, and do not require any special treatment. The computing times within the tables are given in seconds, and are counting •For the iterative method: the time needed for reaching convergence, or, in non-convergent experiments, until a maximum of 2000 iterations have been carried out. •For the direct method: the time spent in defining the matrix plus the time devoted to solving the linear system. When the seed for the iterative method is said to be computed with the (11,10)-direct method, what is meant is that one uses the direct method with I= 11, N= 10 to solve the full problem, and then computes the seed, from the restriction of this solution to Q0, by means of linear interpolation. Lastly, the suffixes “IT” and “D” in the tables stand, respectively, for “iterative” and “direct”. 5.1 Test 1 This is a test taken from [11]. One chooses α(µ, z) = |sin(12µz)|,(26) σ(µ, z) = 1 + sin(12µz) cos(12µz),(27) ψ(µ, z) = ln(2 + µ2+z3),(28) and W,fand gare computed so that ψbe the exact solution. 5.1.1 Results obtained with the iterative method Let us analyze some numerical results obtained with the iterative method. Results collected in Tables 1 and 2 show the influence that the choice of the seed has got in the number of iterations needed for convergence, while Table 3 shows that the number of iterations grows if no relaxation is done. The meaning of “no convergence” in these three Tables 1, 2 and 3 is that convergence has not been reached yet, but the algorithm behaves in a convergent way; in other words, convergence is expected if more iterations are allowed. By way of example, we report in the last row of Table 2 the value of the residual after 2000 iterations. (I, N)Eabs(Q)-IT iterations time-IT (s) (11,10) 7.05 ×10−342 0.17 (33,29) 5.78 ×10−4136 1.83 (101,91) 5.84 ×10−5405 19.01 (0.047 s/it) (321,281) 4.73 ×10−61191 231.68 (0.195 s/it) (1001,901) 2.34 ×10−4No convergence in 2000 iterations. 2212.43 (1.106 s/it) Table 1: Numerical results for test 1 by means of the iterative method with seed given by Equation (14), ω= 2, and ε= 10−8. 9