Elsevier Editorial System(tm) for Applied Mathematics and Computation Manuscript Draft Manuscript Number: Title: Asymptotical Computations for a Model of Flow in Concrete Article Type: Full Length Article Keywords: Flow in concrete; Interface problem; Adaptive mesh selection; Asymptotic theory. Corresponding Author: Dr Othmar Koch, Corresponding Author's Institution: First Author: Pierluigi Amodio, Professor Order of Authors: Pierluigi Amodio, Professor; Chris J Budd, Professor; Othmar Koch; Giuseppina Settanni, Dr.; Ewa B Weinmüller, Professor Abstract: We discuss an initial value problem for an implicit second order ordinary differential equation which arises in models of flow in concrete. Depending on the initial condition, the solution features a sharp interface with derivatives which become numerically unbounded. By using an integrator based on finite difference methods and equipped with adaptive step size selection, it is possible to compute the solution on highly irregular meshes. In this way it is possible to verify and predict asymptotical theory near the interface with remarkable accuracy. Cover Letter Univ.-Doz. Dr. Othmar Koch Institute for Analysis and Scientific Computing (E101) Vienna University of Technology Wiedner Hauptstrasse 8-10 A-1040 Wien AUSTRIA To the Editors of Applied Mathematics and Computation Please find enclosed our manuscript Asymptotical Computations for a Model of Flow in Concrete by P. Amodio, C.J. Budd, O. Koch, G. Settanni, and E. Weinmüller which we would like to submit for possible publication in Applied Mathematics and Computation. Best regards P. Amodio, C.J. Budd, O. Koch, G. Settanni, and E. Weinmüller *Manuscript Click here to view linked References Asymptotical Computations for a Model of Flow in Concrete P. Amodio a, C.J. Budd b , O. Koch c,∗, G. Settanni a, E. Weinmüller c a Dipartimento di Matematica, Università di Bari via E. Orabona 4, 70125 Bari, Italy b Mathematical Sciences, University of Bath Bath BA2 7AY, United Kingdom c Institute for Analysis and Scientific Computing (E101), Vienna University of Technology, Wiedner Hauptstrasse 8–10, A-1040 Wien, Austria Abstract We discuss an initial value problem for an implicit second order ordinary differential equation which arises in models of flow in concrete. Depending on the initial condition, the solution features a sharp interface with derivatives which become numerically unbounded. By using an integrator based on finite difference methods and equipped with adaptive step size selection, it is possible to compute the solution on highly irregular meshes. In this way it is possible to verify and predict asymptotical theory near the interface with remarkable accuracy. Key words: Flow in concrete; Interface problem; Adaptive mesh selection; Asymptotic theory. 1991 MSC: 65L05; 65L12; 34E99. 1 Introduction A model for the time dependent flow of water through a variably saturated porous medium with exponential diffusivity, such as soil, rock or concrete is ∗ Corresponding author. Email addresses: amodio@dm.uniba.it (P. Amodio), cjb@maths.bath.ac.uk (C.J. Budd), othmar@othmar-koch.org (O. Koch), settanni@dm.uniba.it (G. Settanni), e.weinmueller@tuwien.ac.at (E. Weinmüller). Preprint submitted to Elsevier 20 July 2012 given by ! ∂u ∂u ∂ D(u) , (1) = ∂t ∂x ∂x where u(x, t) is the saturation and is the fraction by volume of the pore space occupied by the liquid a distance x into the porous medium, and D(u) = D0 eβu . In this problem the bulk of the liquid resides in the interval x ∈ [0, x∗ (t)] where the moving interface x∗ (t) is called the wetting front, and u ≪ 1 if x > x∗ . The physical derivation of this equation is given in [7], [11], [13]. The numerical treatment of this problem was first discussed in [15]. A comprehensive overview of numerical methods for flow in porous media is given for instance in [10]. In the present paper we do not focus on a simulation of the full model, but adopt a numerical approach to investigate the asymptotical behavior of self-similar solutions of the equation (1), which are stable attractors and take the form u(x, t) = ψ(y), y = x/t1/2 , 0 < y < ∞. If we set θ(y) = eβψ(y) and make a trivial rescaling, it then follows that θ(y) satisfies the ordinary differential equation problem θ(y)θyy (y) = −yθy (y), θy (0) = −γ = −γ(β) < 0, θ(0) = 1, y > 0. (2) The purpose of this paper is to make a numerical study of the solutions of (2) in the limit of large γ which corresponds to a problem with β ≫ 1 with large diffusion when u is not small. A plot of the solution for γ = 3 is given in Figure 1. In this plot we can see that for smaller values of y the solution θ(y) decreases almost linearly, coming close to zero at the point y ≈ 1/γ. As θ approaches zero it follows from (2) that θyy increases. Indeed, we can see from the figure that the solution has high positive curvature close to this point. For larger values of y the solution asymptotes to a constant value so that θ(y) → θ∞ as y → ∞. Of interest are both θ∞ and the point y ∗ at which θyy takes its maximum value, so that θyyy (y ∗) = 0. In analogy to the time-dependent problem we will refer to the latter as the 2 wetting front. Indeed, if y < y ∗ then θ = O(1) and if y > y ∗ then θ = O(e−γ ). An asymptotic study of this problem for large values of γ is given in [8], the 2 results of which are summarised in Section 1. One purpose of this paper is to use an accurate high order numerical method to examine the asymptotic predictions of the latter paper and to accurately determine both the location of the wetting front and the value of θ∞ . This problem poses a significant numerical challenge. Whilst the error terms in the asymptotic expansions (say of y ∗ ) decay polynomially in γ as γ increases, it is shown in [8] that θyy (y ∗ ) ∼ eγ 2 2 and θ∞ ∼ e−γ . Hence, for even moderately large values of γ we see that the solution features a near-singularity in its higher derivatives, and moreover the solution is almost identically equal to 0 for larger values of y. Indeed in the calculations reported in this paper we have values of θyy ∼ 10138 and θ ∼ 10−141 . Conventional initial value solvers such as the Gear solvers in the Matlab routine ode15s are unable to deal with the form of these singularities for values of γ at which the asymptotic results are expected to become sharp (for example values of γ > 10). For the larger values of γ required to study the asymptotic results a different numerical approach is required, and it is the difficulties described above that put this problem in the focus of the higher order methods that we will discuss in this paper. The numerical challenge is thus to find discretization schemes and adaptive mesh selection methods which are robust with respect to the dramatic variation in the step sizes which may be expected due to the predicted solution behavior as the curvature of θ rapidly increases. Accordingly, in this paper we propose the construction of a variable step size numerical scheme directly for the second order problem (2) without resorting to a transformation to first order form. The idea of discretizing each derivative of the second order continuous problem by means of finite difference schemes has been largely considered in the past, in particular in the numerical solution of partial differential equations using schemes of low order (and with small stencils). In [4] and later in [1,2,5] it was proposed to apply high order finite difference schemes for the solution of boundary value problems for ordinary differential equations by considering different formulae with the same order in the initial and final points of the grid, following the idea inherited by boundary value methods [6]. One major advantage of this approach is associated with the fact that the vector of unknowns contains only the solution of the problem. This choice, on the one hand, reduces the computational cost of the algorithm with respect to the standard approaches which transform the second order equation into a first order one containing both the solution and its first derivative as unknowns. On the other hand, it simplifies the stepsize variation and allows the solution of fully implicit problems without any change in the approach. Moreover, in [12] it was demonstrated that a solution in the second 3 order formulation may yield an advantage with respect to the conditioning of the linear systems of equations associated with the discretization scheme. Also, it is clear that in this approach, only the smoothness of the solution, but not of its derivative, affects the step-size selection process, see [9]. Finally, in the original formulation there is complete freedom in the choice of the schemes approximating each derivative which could not only depend on whether an initial or boundary value problem is solved, but also on the problem data and the discretization points. In [3] the same idea is applied to IVPs for ODEs. Two approaches are proposed to take care of the first derivative in the initial point: we can choose to approximate this initial condition by means of an appropriate formula or to define difference schemes which make use of the analytical first derivative in the first point. As a general belief, the second approach seems to be preferable but, for the flow in concrete problem the first strategy is used, since the solution’s derivative may be quite different from the solution itself (and hence affect the stepsize variation). The computations which enable to provide the mentioned results are only possible due to this highly flexible, adaptive finite difference method which is able to cope with extremely unsmooth step-size sequences. The present numerical investigation is successful. In particular for the case where γ ≫ 1, the simulations give very stable results which closely match the theoretical results and the expected errors, both in confirming certain of the asymptotic predictions for y ∗ and θ∞ derived in [8] and also in giving strong evidence for further asymptotic results which are beyond the existing theory. This lends validity to both the numerical method and the asymptotic calculation. The layout of the remainder of this paper is as follows. In Section 2 of the paper, we briefly describe the asymptotic theory for (2) given in [8] and compare the results with our numerical simulations for values of γ ≈ 18 for which θyy (y ∗ ) ≈ 10138 . We show not only that the numerical results give good agreement with the existing results, but moreover the computations enable us to determine constants that are not provided by the a priori analysis. The numerical method itself is described in Section 3. In Section 4 we show the form of the solution θ(y) and give plots of the solution and of its derivatives to illustrate the origin of the numerical difficulties described above. In Section 5 a more detailed analysis of the solution behaviour is given in the so-called mid-range of the independent variable in which y ≈ y ∗ and θyy (y) has its very high values. Finally we draw some conclusions from this work. 4 2 The asymptotic properties of the solution of (2) for large γ. In this section we outline the main asymptotic predictions for the solutions of (2) which are presented in [8] and will show how these are supported by the numerical calculations. In [8] a matched asymptotic expansion method is used to find an asymptotic form of both the wetting front y ∗ and the final value of log(θ∞ ), expressing both as asymptotic series developed in powers of 1/γ with γ ≫ 1. The asymptotic series are found by matching descriptions of the solution in an inner range with y < y ∗ ≈ 1/γ, a mid-range with y ≈ y ∗ and an outer range with y ≫ y ∗ . It is in the mid-range where the most delicate behaviour occurs, with exponentially large (in γ) values for θyy making the asymptotic theory delicate in this case. We will describe the precise form of the solution in the mid-range in Section 4. To summarize: In [8] it is proposed that as γ → ∞ the value of θ(1/γ) is given asymptotically by the expression 1 θ γ ! = 1 log(γ) b − + 4, 2 4 2γ γ γ (3) where 1 b=b +O γ2 ∗ ! with b∗ = 11 1 − log(2) ≈ 0.57009307.... 12 2 (4) Similarly, the final value θ∞ satisfies the asymptotic relation log(θ∞ ) = −γ 2 − α 1 + 2, 2 γ (5) where it is conjectured that α takes the form ! 1 α =α +O , γ2 ∗ (6) with α∗ an unknown constant. Furthermore, the location of the wetting front y ∗ is given asymptotically by the expression y∗ = where 1 β 1 + 3 + 5, γ 2γ γ ! 1 , β =β +O γ2 At the wetting front we have ∗ ∗ θ(y ) ≈ eθ∞ with β ∗ = ∗ and θyy (y ) ≡ θ̂yy ≈ 5 (7) 11 . 12 eγ (8) 2 −1/2 γ2 . (9) Using the finite difference schemes described in the next section, we study these asymptotic results by solving the initial value problem, θθyy = −yθy , θy (0) = −γ, θ(0) = 1, (10) for γ = 2, 3, . . . 18. For γ > 18 the solution for large y in which θ ≈ θ∞ = 2 10−141 ∼ e−γ is too small to be computed accurately and for γ > 26 it is smaller than the smallest positive double precision machine number. Note that although γ = 18 poses a serious challenge for the numerical method, the predicted asymptotic error of order 1/γ 2 is not unreasonably small. We consider methods of orders p equal to 4, 6, 8, and 10 to understand how the order of the numerical method influences the accuracy of the approximate solutions. For the purpose of illustration, we will generally give the results for p = 4 and p = 8. As a general remark, we observe that the discretization errors are larger than the round-off errors of the floating point operations. Moreover, we propagate the solution using variable stepsizes, and for a fixed tolerance, methods of higher order require smaller grid size but they do not necessarily achieve better precision. This is most probably due to the fact that higher derivatives of the solution θ are extremely unsmooth. The solution has been computed on finite intervals [0, y∞ ], with the interval endpoint y∞ = 10, 20, 30, 40, 50, to see how strongly the value of θ(y∞ ) depends on the length of the interval of integration. It turns out that, especially for large values of γ, the influence of the interval length is negligible. Therefore, from now on, we use the interval [0, 10] and θ(∞) ≈ θ∞ := θ(10) for all calculations. In Table 1 we present the results of the numerical computations of the various terms in the expressions above. All data is given for computations with γ = 10 and γ = 18, respectively. Since the reference values are asymptotically correct for large values of γ, we expect to observe higher accuracy for larger γ. In Rows 1 and 2 of Table 1, we specify the values of θ∞ and y ∗, where y ∗ is defined to be the point where the numerically computed value of θyy takes its maximal value in the interval [0, y∞ ] = [0, 10]. In Row 3 we first report the relative error in approximating − log(θ∞ ) by γ 2 + 21 and then in Row 4 we compute an approximation for α as defined in (5). The values in Row 3 confirm that the value of θ∞ has been computed accurately (the relative error is always smaller than γ −4 ) and the computed value for α in expression (5) is α ≈ −0.08. We define ! 1 1 , s =γ y −γ = +O 2 γ2 ∗ 3 ∗ 2 (11) and in row 5 we specify the relative error |s∗ − 0.5|/0.5 and then in Row 6 we use (7) to approximate the value of β. Again the relative error is smaller than γ −2 while the value of β is close to 0.95. 6 Rows 7 and 8 contain the values of θ(y ∗ ) and θyy (y ∗), showing the very rapid increase in the value of the latter with γ. In Rows 9 and 10 we specify the relative errors, θyy (y ∗ ) − θ̂yy |θ(y ∗ ) − eθ∞ | , (12) , eθ∞ θ̂yy respectively. Finally, Rows 11 and 12 contain the numerical value of b as given in (3) and the relative error. In all of these expressions the relative error decreases as γ increases, demonstrating that to leading order both the numerical and asymptotic calculations are correct. Given the relatively modest values of γ used in the computations, and the corresponding large asymptotic errors, it is necessary to post-process the numerical results to determine the finer structure of the solution and to verify the accuracy of the methods used. This post-processing allows us to make some further conjectures about the asymptotic solution. Taking the asymptotic expressions for α, β and b in expressions (5), (6) and (3) in the precise form A , γ2 B β = β∗ + 2 , γ C b = b∗ + 2 , γ α = α∗ + (13) (14) (15) we substitute the results from Table 1 for γ = 10 and γ = 18 to find the corresponding values of α∗ , β ∗ , b∗ , A, B, C. This calculation gives α∗ ≈ −0.0834 ≈ −1.0004/12, with ∗ β ≈ 0.9163 ≈ 10.9954/12, with ∗ b ≈ 0.5620 with C ≈ 4.81. A ≈ −0.089, B ≈ 2.96, These results are close to the asymptotical results given in (8), (4) for which β ∗ = 11/12 ≈ 0.9166666, and b∗ ≈ 0.57009, while the result for α∗ strongly implies the new result that 1 α∗ = − . 12 We conclude from this numerical calculation that there is strong evidence for the correctness of the asymptotic expressions in [8] and that the wetting front y ∗ is located at the point y∗ ≈ 1 11 2.96 1 + 3+ + 7 . 5 γ 2γ 12γ γ 7 Moreover, 1 1 0.089 − − . 2 2 12γ γ4 log(θ∞ ) ≈ −γ 2 − 3 The Numerical Method Here we describe the finite difference schemes underlying the computations presented in this manuscript. These have proven robust with respect to the unusually strong variations in the stepsizes encountered close to the points of high curvature experienced in this problem. The methods are applied to second-order initial value problems in their original formulation, and use different formulae for the appearing derivatives. The difference formulae are defined for equidistant grids and used on small subintervals of the time integration. Stepsize variation is performed after each subinterval as in the general block– BVM framework [6]. For the sake of clarity we consider the numerical solution of a general second order initial value problem  f (y, θ, θ ′ , θ ′′ ) θ(y 0) = θ0 , = 0, θ′ (y0 ) = θ0′ . (16) We first specify an initial stepsize h0 and a grid of equispaced points Y = [y0 , y1 , . . . , yn ], yi = y0 + ih0 . Denote the corresponding numerical approximation by Θ = [θ0 , θ1 , . . . , θn ]. Following the idea in [3–5], we discretize separately the derivatives in (16) by means of suitable high order finite difference schemes θ (ν) (yi ) ≃ (ν) θi r 1 X (s,ν) = ν αs+j θi+j , h j=−s (s,ν) (17) where ν = 1, 2 represents the derivative index and αs+j are the coefficients of the method which are fixed in order to reach the maximally possible order of accuracy. The integers s and r represent the number of left and right values required to approximate θ(ν) (yi ), and are strictly related to the order and the stability of the formula. For this problem we choose, when possible, r = s obtaining formulae (called ECDF in [5]) of even order p = 2s for both the first 8 Table 1 Asymptotic results: Part 1. γ = 10 Order 4 8 θ∞ 2.254440321030590e−44 2.254439903035485e−44 y∗ 1.005094588206025e−01 1.005094588173695e−01 | log(θ∞ ) + γ 2 + 0.5| γ 2 + 0.5 8.381534862525284e−06 8.383379735339227e−06 α −8.423442536837911e−02 −8.425296634015922e−02 |s∗ − 0.5|/0.5 1.891764120503581e−02 1.891763473901165e−02 β 9.458820602506902e−01 9.458817369495821e−01 θ(y ∗ ) 6.126853837541197e−44 6.179228381160976e−44 θyy (y ∗ ) 1.648468426869335e+41 1.648412414631253e+41 2.203452148159067e−04 8.326316736789704e−03 1.106641417450169e−02 1.103205980544998e−02 b 5.138760976981811e−01 5.138757975462900e−01 |b − b̂|/b̂ 9.861017615723677e−02 9.861070265352932e−02 |θ(y ∗ ) − eθ∞ | eθ∞ |θyy (y ∗ ) − θ̂yy |/θ̂yy γ = 18 Order 4 8 θ∞ 1.178498689884995e−141 1.178497239286140e−141 5.564177918914028e−02 5.564177918834908e−02 | log(θ∞ ) + + 0.5| γ 2 + 0.5 7.913110850624213e−07 7.951042681243711e−07 α −8.319686486129285e−02 −8.359567254206013e−02 |s∗ − 0.5|/0.5 5.712462132123619e−03 5.712452903708254e−03 β 9.254188654173835e−01 9.254173704034648e−01 θ(y ∗ ) 3.198906899145607e−141 3.211721064871286e−141 θyy (y ∗ ) 9.664468728333972e+137 9.664458760222388e+137 1.431149209001307e−03 2.570147096053791e−03 3.363040886020695e−03 3.362005998844594e−03 b 5.471501619326116e−01 5.471487116780027e−01 |b − b̂|/b̂ 4.024415556753806e−02 4.024669945847258e−02 y∗ γ2 |θ(y ∗ ) − eθ∞ | eθ∞ |θyy (y ∗ ) − θ̂yy |/θ̂yy 9 and the second derivative. For example, we have for Order 4: 1 h2 θ′′ (yi ) ≈ − 12 θi−2 + 34 θi−1 − 52 θi + 34 θi+1 − h θ′ (yi ) ≈ 1 θ 12 i−2 − 32 θi−1 + 23 θi+1 − 1 θ , 12 i+2 1 θ , 12 i+2 Order 6: h2 θ′′ (yi ) ≈ 1 θ 90 i−3 − 3 θ 20 i−2 1 h θ′ (yi ) ≈ − 60 θi−3 + + 23 θi−1 − 3 θ 20 i−2 49 θ 18 i + 23 θi+1 − − 34 θi−1 + 43 θi+1 − 3 θ 20 i+2 3 θ 20 i+2 + + 1 θ , 90 i+3 1 θ . 60 i+3 We observe that the coefficients are symmetric and skew-symmetric for the second and first derivatives, respectively. We have successfully used these schemes up to the order 10. We compute the unknowns in Θ by solving the nonlinear system (1) (2) f (yi , θi , θi , θi ) = 0, i = 1, . . . , n − 1, (18) together with the initial conditions. Since the above symmetric formulae of order p > 2 cannot be used to approximate θ(ν) (yi ), i = 1, . . . , p/2 − 1 and i = n − p/2 + 1, . . . , n − 1, schemes in (17) with different stencils (but the same order) must be provided at the beginning and the end of the grid. We call them initial and final formulae (see [6]). Examples of the initial schemes 10 are Order 4 5 5 1 7 1 1 h2 θ′′ (y1 ) ≈ θ0 − θ1 − θ2 + θ3 − θ4 + θ5 , 6 4 3 6 2 12 1 5 3 1 1 h θ′ (y1 ) ≈ − θ0 − θ1 + θ2 − θ3 + θ4 . 4 6 2 2 12 Order 6 7 7 27 19 67 9 1 11 θ0 − θ1 − θ2 + θ3 − θ4 + θ5 − θ6 + θ7 , 10 18 10 4 18 5 2 180 107 21 13 17 3 4 1 11 θ0 + θ1 − θ2 + θ3 + θ4 − θ5 + θ6 − θ7 , h2 θ′′ (y2 ) ≈ − 180 90 10 18 36 10 45 90 1 77 5 5 5 1 1 h θ′ (y1 ) ≈ − θ0 − θ1 + θ2 − θ3 + θ4 − θ5 + θ6 , 6 60 2 3 6 4 30 1 2 7 4 1 2 1 h θ′ (y2 ) ≈ θ0 − θ1 − θ2 + θ3 − θ4 + θ5 − θ6 . 30 5 12 3 2 15 60 h2 θ′′ (y1 ) ≈ The final schemes used to approximate θ(ν) (yn ) and, for the order 6, θ(ν) (yn−1 ), are not reported. Anyway, the coefficients of these methods correspond to the initial ones given above, but in reversed order and with the opposite sign in case of the first derivative. Note that for the second derivative, the order of the initial methods is p = r + s − 1 (we need one additional value to obtain the required order). The number of grid points n ≥ p computed by solving (18) is linked to stability and computational cost. Larger n means better stability properties but higher computational cost [3,6]. For this problem we have fixed n = p + 4. The structure of the coefficient matrix associated with the second derivative for p = 10 is depicted in Figure 2. We illustrate in this example 4 initial and 4 final schemes, and 5 times the main scheme. The size of the matrix is (n−1)×(n+1), but the first column is multiplied by the known value θ0 . θ′ (0) is approximated by a suitable starting scheme obtained by choosing s = 0 in (17). Examples of these last formulae are: Order 4: 25 hθ′ (y0 ) ≈ − 12 θ0 + 4θ1 − 3θ2 + 34 θ3 − 41 θ4 = hθ0′ , Order 6: θ + 6θ1 − hθ′ (y0 ) ≈ − 49 20 0 15 θ 2 2 + 11 20 θ 3 3 − 15 θ 4 4 + 65 θ5 − 16 θ6 = hθ0′ . Once the solution in yi = y0 + ih0 , i = 1, . . . , n, has been approximated, the stepsize is changed according to an estimate of the local error based on the mesh halving strategy and the algorithm is iterated on a subsequent subinteval. The code solving the flow in concrete problem uses a classical time stepping strategy [14] with safety factor 0.7. Since however the second last values θn−1 ′ and θn−1 are better approximated than θn and θn′ , we discard the last value in ′ Θ, use formula (17) with i = n − 1, r = 1 to compute θn−1 , and continue the ′ algorithm with the initial conditions [θn−1 , θn−1 ]. In order to obtain the results in this paper, for all the even orders from 4 to 10 we have used initial stepsize h0 = 1e − 3 and error tolerance tol = 1e − 12 for both the solution and the Newton iteration. 4 Solution plots In Fig. 3, we show the graphs in a logarithmic scale of the numerical solution and its first and second derivatives for the value γ = 18 computed by the numerical method of order 8. The rapid change in θ and its derivatives close to the wetting front is very clear from these figures. It was found that the approximation of the first and the second derivative becomes unreliable in the area where the solution becomes constant. In fact, since derivatives are computed by means of linear combinations of the solution values (using finite difference formulae (17)), it is really unlikely that they can be lower (in absolute value) than θ∞ EPS, where EPS is the machine precision, and in this case must be treated as 0. Secondly, in Fig. 4 we plot the variation of the step sizes for the methods of orders 4, 6, 8, and 10. The used tolerance is 1e−12 and the initial step size is 1e−3. We note that the orders 6 and 8 require the smallest number of meshpoints, since the size of each block is p + 4. 1 We observe the very small step sizes used in these methods close to the point of high curvature. Note that the same minimal stepsize was defined for all orders, and this stepsize is reached with fewer but larger steps if the order of the method is higher. 5 Asymptotical properties revisited: Mid–range calculation Next, we consider the asymptotical solution properties in the regime where y lies in the mid-range close to y ∗ ≈ 1/γ. It is convenient to rescale both the 1 Clearly, the number of meshpoints is equal to p + 4 times the number of blocks (see Section 3). 12 dependent and independent variables in this region according to (11), such that 1 s y = + 3 and v(s) = γ 4 θ(y), γ γ with |s| varying between 1 and γ 2 . We first discuss the case of negative s. Case 1: Let −γ 2 ≪ s ≪ −1. In [8] an asymptotic series for the function v(s) is developed in the form ! 1 . v(s) = γ v0 (s) + v1 (s) + O γ2 2 (19) A careful application of the method of matched asymptotic expansions then implies that 1 (20) v0 (s) = − s 2 and 1 1 log(2) v1 (s) = 2 log(γ)s − log(γ) + b − s − log −s + , 2 2 2  where b =    (21) 11 1 − log(2). 12 2 In Figures 5 and 6, we plot the values of v1 (s) and their numerical approximations together with the relative errors (logarithmic scale), respectively. Computations were realized with the method of order 8. We can see the desired behaviour in the asymptotical regime −γ 2 ≪ s < 0. Note furthermore that as expected the approximation improves as γ increases. We next consider the mid-range with positive values of s. It is in this range that we see the rapid transition from polynomial to exponential decay. Case 2: Let 1 ≪ s ≪ γ 2 . The asymptotic form of v(s) is given by [8],  v(s) = v∞ 1 + e−γe e−(s−s ∗ )/v ∞  , (22) where s∗ = 12 + O(1/γ 2) and γe ≈ 0.577215665 is the Euler–Mascheroni constant. For a large value of γ the value of v∞ is very small, and therefore, v(s) in (22) is constant and coincides with v∞ ≡ θ∞ /γ 4 . Figures 7 and 8 show the graphs of v(s) and their absolute errors obtained based on (22), computed by the method of order 8 for γ = 10 and γ = 18, 13 respectively. Note that in spite of the small solution values, the numerical results are still meaningful with an accuracy of about five digits. However, the very small value of the solution makes all comparisons difficult. Conclusions We have discussed the numerical solution of a second–order ODE problem which arises as a model of water flow through concrete. The solution of this problem features an interface with very high values of the higher derivatives. An adaptive finite difference method is employed to approximate the solution numerically and verify predictions derived by asymptotic theory. The solution algorithm can cope with a large variation in the step-size and thus can serve to accurately approximate the location of the interface and asymptotic results about the solution structure. Thus the asymptotic results are verified, parameters not given by the a priori analysis are determined and new predictions about the solution structure are indicated. References [1] P. Amodio and G. Settanni. A deferred correction approach to the solution of singularly perturbed BVPs by high order upwind methods: implementation details. AIP Conf. Proc., 1168:711–714, 2009. [2] P. Amodio and G. Settanni. Variable step/order generalized upwind methods for the numerical solution of second order singular perturbation problems. JNAIAM J. Numer. Anal. Indust. Appl. Math., 4:65–76, 2009. [3] P. Amodio and G. Settanni. High order finite difference schemes for the solution of second order initial value problems. JNAIAM J. Numer. Anal. Indust. Appl. Math., 5:3–16, 2010. [4] P. Amodio and I. Sgura. High-order finite difference schemes for the solution of second-order BVPs. J. Comput. Appl. Math., 176:59–76, 2005. [5] P. Amodio and I. Sgura. High order generalized upwind schemes and numerical solution of singular perturbation problems. BIT, 47:241–257, 2007. [6] L. Brugnano and D. Trigiante. Solving Differential Problems by Multistep Initial and Boundary Value Methods. Gordon and Breach Science Publishers, Amsterdam, 1998. [7] W. Brutsaert. Universal constants for scaling the exponential soil water diffusivity. Water Resour. Res., 15(2):481–483, 1979. 14 [8] C.J. Budd and J. Stockie. Asymptotic behaviour of wetting fronts in porous media with exponential moisture diffusivity. University of bath report, University of Bath, 2012. [9] J. Cash, G. Kitzhofer, O. Koch, G. Moore, and E.B. Weinmüller. Numerical solution of singular two-point BVPs. JNAIAM J. Numer. Anal. Indust. Appl. Math., 4:129–149, 2009. [10] Z. Chen, G. Huan, and Y. Ma. Computational Methods for Multiphase Flows in Porous Media. SIAM, Philadelphia, PA, 2006. [11] B. E. Clothier and I. White. Measurement of sorptivity and soil water diffusivity in the field. Soil Sci. Soc. Amer. J., 45:241–245, 1981. [12] G. Kitzhofer, O. Koch, P. Lima, and E. Weinmüller. Efficient numerical solution of the density profile equation in hydrodynamics. J. Sci. Comput., 32:411–424, 2007. [13] C. Leech, D. Lockington, and P. Dux. Unsaturated diffusivity functions for concrete derived from NMR images. Mater. Constr., 34:413–418, 2003. [14] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. Numerical Recipes in C — The Art of Scientific Computing. Cambridge University Press, Cambridge, U.K., 1988. [15] M. Rose. Numerical methods for flows through porous media I. Math. Comp., 40(162):435–467, 1983. 15 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Fig. 1. The solution for γ = 3. In this figure we can see the initial linear decrease of the solution towards zero, the point of high curvature for y = y ∗ ≈ 1/3 and the very small value assumed for y > y ∗ . Fig. 2. Structure of the coefficient matrix approximating the second derivative, p = 10 and n = p + 4. 16 0 10 −50 θ(y) 10 −100 10 −150 10 −4 10 −3 10 −2 −1 10 10 0 10 1 10 y (a) Numerical solution. 50 10 0 10 −50 θ’(y) 10 −100 10 −150 10 −200 10 −4 10 −3 10 −2 −1 10 10 0 10 1 10 y (b) Numerical approximation of minus the first derivative. 200 10 100 θ’’(y) 10 0 10 −100 10 −200 10 −4 10 −3 10 −2 −1 10 10 0 10 1 10 y (c) Numerical approximation of the second derivative. Fig. 3. Solution and its derivatives for γ = 18. Missing points correspond to negative values of −θy (y) and θyy (y) 17 0 0 10 10 −50 −50 Step size 10 Step size 10 −100 −100 10 10 −150 10 −150 0 1000 2000 3000 4000 5000 Number of blocks 6000 7000 10 8000 0 500 (a) order = 4. 1000 1500 Number of blocks 2000 2500 (b) order = 6. 0 0 10 10 −50 −50 Step size 10 Step size 10 −100 −100 10 10 −150 −150 0 200 400 600 800 1000 Number of blocks 1200 1400 10 1600 0 200 (c) order = 8. 400 600 800 1000 1200 Number of blocks 1400 1600 (d) order = 10. Fig. 4. Step size variation for each block, γ = 18. 2 10 Numerical Theoretical 1 Numerical Theoretical 0 0 −1 −10 v (s) −3 1 v1(s) −2 −20 −4 −30 −5 −6 −40 −7 −8 −16 −14 −12 −10 −8 s −6 −4 −2 −50 −100 0 −80 −60 −40 −20 0 s (a) γ = 4. (b) γ = 10. 10 20 Numerical Theoretical 0 Numerical Theoretical 0 −20 −10 −40 −20 −60 1 v (s) −30 1 v (s) 10 −40 −80 −100 −50 −120 −60 −140 −70 −80 −150 −160 −100 −50 −180 −350 0 s −300 −250 −200 −150 −100 −50 0 s (c) γ = 12. (d) γ = 18. Fig. 5. Analytical values and numerical approximation for v1 , order 8. 18 1800 3 3 10 10 2 10 2 10 1 Error Error 10 1 10 0 10 0 10 −1 10 −1 10 −16 −2 −14 −12 −10 −8 s −6 −4 −2 10 −100 0 −80 −60 −40 −20 0 s (a) γ = 4. (b) γ = 10. 3 3 10 10 2 10 2 10 1 10 1 Error Error 10 0 10 0 10 −1 10 −1 10 −2 10 −2 −3 10 −150 −100 −50 10 −350 0 −300 −250 −200 s −150 −100 −50 0 s (c) γ = 12. (d) γ = 18. Fig. 6. Relative error of v1 , order 8. γ = 10 order = 8 −40 2.2545 x 10 γ = 10 order = 8 −45 2 x 10 2.2545 1.5 v(s) Error v(s) 2.2545 2.2545 0.5 2.2545 2.2545 1 0 20 40 60 80 0 100 s 0 20 40 60 80 s (a) v(s), max |v(s)| = 2.254510924061732e − 040 s Fig. 7. γ = 10, order 8 19 (b) Absolute error of v(s) 100 γ = 18 order = 8 −136 1.2373 x 10 3 1.2373 2.5 1.2373 2 v(s) Error v(s) 1.2373 1.2373 1 1.2373 0.5 0 50 100 150 200 250 300 0 350 x 10 1.5 1.2373 1.2373 γ = 18 order = 8 −141 3.5 0 s 50 100 150 200 250 300 s (a) v(s), max |v(s)| = 1.237327181687361e − 136 s Fig. 8. γ = 18, order 8 20 (b) Absolute error of v(s) 350