A priori estimation of a time step for numerically solving parabolic problems

This work deals with the problem of choosing a time step for the numerical solution of boundary value problems for parabolic equations. The problem solution is derived using the fully implicit scheme, whereas a time step is selected via explicit calculations. The selection strategy consists of the following steps. First, using the explicit scheme, we calculate the solution at a new time level. Next, we employ this solution in order to obtain the solution at the previous time level (the implicit scheme, explicit calculations). This solution should be close to the solution of our problem at this time level with a prescribed accuracy. Such an algorithm leads to explicit formulas for the calculation of the time step and takes into account both the dynamics of the problem solution and changes in coefficients of the equation and in its right-hand side. The same formulas for the evaluation of the time step we get using a comparison of two approximate solutions, which are obtained using the explicit scheme with the primary time step and the step that is reduced by half. Numerical results are presented for a model parabolic boundary value problem, which demonstrate the robustness of the developed algorithm for the time step selection.


Introduction
In numerically solving boundary value problems for time-dependent equations, emphasis is on discretizations in time [1,2,7]. For parabolic equations of second order, unconditionally stable schemes are constructed using implicit approximations [9,10,11]. Two-level schemes are commonly used in computational practice, whereas multilevel schemes occur more rarely. For unconditionally stable schemes, a time step is selected only due to the accuracy of the approximate solution.
The problem of the control over a time step is relatively well resolved for the numerically solving Cauchy problem for systems of differential equations [3,5,6]. The basic approach involves the following stages. First, we perform additional calculations in order to estimate the error of the the approximate solution at a new time level. Further, a time step is estimated using the theoretical asymptotic dependence of accuracy on a time step. After that we decide is it necessary to correct the time step and to repeat calculations.
Additional calculations for estimating the error of the approximate solution may be performed in a different way. In particular, it is possible to obtain an approximate solution using two different schemes that have the same theoretical order of accuracy. The most famous example of this strategy involves the solution of the problem on a separate time interval using a preliminary step (the first solution ) and the step reduced by half (the second solution). In numerically solving the Cauchy problem for systems of ordinary differential equations, there are are also applied nested methods, where two approximate solutions of different orders of accuracy are compared.
In the above-mentioned methods of selecting a time step, a posteriori estimation of accuracy is employed. In this case, we decide is this time step acceptable or it is necessary to change it for re-calculations (increase or reduced and how much) only after performing calculations. Such strategies can be also applied to the approximate solution of unsteady boundary value problems using a more advanced a posteriori analysis [4,8,13].
In this paper, we consider an a priori selection of a time step for the approximate solution of boundary value problems for parabolic equations. To obtain the solution at a new time level, the backward Euler scheme is employed. The time step at the new time level is explicitly calculated using two previous time levels and takes into account changes in the equation coefficients and its right-hand side. The paper is organized as follows. In Section 2, we consider a Cauchy problem for a system of linear ordinary differential equations that is obtained from numerically solving boundary value problems for parabolic equations after discretization in space. For the approximate solution, estimates for stability are presented along with estimates for accuracy in the corresponding Hilbert space. Formulas for the selection of a time step are obtained in Section 3 using a comparison of the problem solutions corresponding to the forward time level and backward one. In Section 4, we show that similar expressions for a time step can be obtained via making a comparison of the solutions derived with one time step and two half steps. Section 5 presents numerical results for a model boundary value problem for a one-dimensional parabolic equation obtained on the basis of the developed algorithm for selecting a time step.

Model problem
Let us consider the Cauchy problem for the linear equation supplemented with the initial condition The problem is investigated in a finite-dimensional Hilbert space H. Assume that in H. Due to the non-negativity of the operator A, for the problem (1), (2), we have the following estimate for stability with respect for the initial data and the right-hand side: The problem (1), (2) results from finite difference, finite volume or finite element approximations (lumped masses scheme [12]) for numerically solving boundary value problems for a parabolic equation of second order. In this problem, an unknown function u(x x x,t) satisfies the equation To solve numerically this time-dependent problem, we introduce a non-uniform grid in time: We will employ notation f n = f (t n ). For the problem (1), (2), we apply the fully implicit scheme, where the transition from the current time level to the next one is performed as follows: starting from the initial condition Under the restriction A n+1 ≥ 0, from (4), it follows immediately that the approximate solution satisfies the level-wise estimate Thus, we obtain the discrete analog of the estimate (3): corresponding to the problem (4), (5). For the error z n = y n − u n of the approximate solution, we have the problem Here ψ n+1 stands for the truncation error: Similarly to (6), we get the estimate for error: Therefore, to control the error, we can employ the summarized error τ n+1 δ over the interval t n ≤ t ≤ t n+1 . In this case, a value δ defines the same level of the error over the entire interval of integration.

Algorithm for estimation of a time step
If we will be able to calculate the truncation error ψ n+1 , then it will be possible to get a posteriori estimate of the error. Comparing ψ n+1 with the prescribed error level δ , this makes possible to evaluate the quality of the choice of the time step τ n+1 . Namely, if ψ n+1 is much larger (smaller) than δ , then the time step is taken too large (small), and if ψ n+1 is close to δ , then this time step is optimal. Thus, we have The problem is that we cannot evaluate the truncation error, since it is determined using the exact solution that is unknown. Because of this, we must focus on some estimates for the truncation error that guarantee the fulfilment of (8).
The following strategy is proposed for the correction of the time step. The step τ n+1 is selected from the conditions: Forward step. Using the explicit scheme, we calculate the solution v n+1 at the time level t n+1 ; Backward step. From the obtained v n+1 , applying the implicit scheme, we determine v n at the time level t n (explicit calculations); Step selecting. The step τ n+1 is evaluated via closeness between v n and y n .
In fact, we carry out the back analysis of the error of the approximate solution over the interval t n ≤ t ≤ t n+1 using two schemes (explicit and implicit) of the same accuracy. Let us present the formulas for selecting a time step. The solution v n+1 is determined from the equation v n+1 − y n τ n+1 + A n y n = f n .
For v n , we have v n+1 − v n From (9), (10), we immediately get The first two terms are associated with the time derivative applied to the problem operator and to the right-hand side. To evaluate them approximately, it seems reasonable to use the time step from the previous time level. But this may be inconvenient to implement.
For instance, we have and therefore we have to evaluate the difference derivative of the right-hand side for t n ≤ t ≤ t n+1 . The problem is that the derivation of such estimates involves the unknown value t n+1 . The simplest approach is to evaluate this derivative using the previous time step: But in this case, if τ n+1 > τ n , then we cannot detect significant changes in the right-hand side for t n + τ n ≤ t ≤ t n + τ n+1 .
To resolve the problem, it is possible to use the standard methods available to control a time step for numerically solving time-dependent problems. The first method restricts the growth of the time step with respect to the previous value. We set where γ > 1 is a numerical parameter. The second requirement is that the step cannot be too small: where τ 0 is a specified minimum time step. Under the assumption (12), we can estimate the time derivative of the right-hand side, putting Therefore For the last term in the right-hand side of (11), in view of (4), we have τ n A n+1 (y n − y n−1 ).
With this in mind, the equality (11) is replaced by the approximate equality: The value of v n − y n we associate with the solution error over the interval t n ≤ t ≤ t n+1 . Because of this, we set v n − y n ≤ τ n+1 δ .
From (15), we have v n − y n τ 2 In view of (12), (13), (16), from (15), we obtain the following formula for calculating the time step: This formula for selecting a time step reflects clearly (see the denominator in the expression for γ n+1 ) corrective actions, which are related to the time-dependence of the problem operator (the first part) and the right-hand side (the second part) as well as to the time-variation of the solution itself (the third part).

Estimation of a time step on the basis of step doubling
To solve numerically the Cauchy problem, the traditional strategy is to select an integration step using a comparison of the approximate solution obtained by the preliminary step with the solution calculated with the step reduced by half. For numerically solving problem (1), (2), we use fully implicit scheme (4), (5). We employ the explicit scheme over the interval t n ≤ t ≤ t n+1 in order to select the time step τ n+1 . The selection strategy includes: Calculation with an integer step. Using the explicit scheme, we determine the solution v n+1 at the time level t n+1 via the step τ n+1 ; Calculation with a half-integer step. Using the explicit scheme, we calculate the solution w n+1 at the time level t n+1 employing the step 0.5τ n+1 ; Step selecting. The step τ n+1 is evaluated through the closeness between v n+1 and w n+1 .

Thus, we arrive at
The right-hand side of (22) coincides with the right-hade side of (13) with an accuracy of a factor. Similarly to (18), we can formulate the rule for selecting the time step: In fact, we have come to the same rule for the estimation of the time step -the factor 4 has not any essential matter.

Numerical experiments
To demonstrate the performance of the proposed algorithm (12), (16) for selecting a time step based on the implicit scheme for solving the problem (1), (2), let us consider the boundary value problem for a one-dimensional parabolic equation. Let u(x,t) satisfies the equation as well as the boundary and the initial conditions: u(x, 0) = u 0 (x), To solve approximately the problem (24)-(26), we apply finite difference discretization in space. Let us introduce a uniform grid with a step h: and ω is the set of interior grid points, whereas ∂ ω is the set of boundary points (ω = ω ∪ ∂ ω). On the set of grid functions such that u(x) = 0, x / ∈ ω, we introduce a Hilbert space H, where the inner product and the norm are defined as: The grid operator A(t) is written as follows: The operator A(t) is self-adjoint, and if p(t) ≥ 0, then it is positive definite in H. Thus, after the spatial discretization of (24)-(26), we arrive to the problem (1), (2). As a test problem, we consider the problem (24)-(26) with T = 0.1 and the discontinuous coefficient p(t) and the discontinuous f (t) defined as follows: The problem is solved on the grid with M = 100, the calculations are performed using the sufficiently small time step τ 1 = τ 0 = 1 · 10 −6 . First, we Consider the case, where the initial condition (26) is taken in the following form: If we specify the error level δ = 0.1 and the parameter γ = 1.5, then the time step produced by the algorithm (14), (18) has the form depicted in Fig. 1. The total number of time steps is N = 268. In this figure, we observe essential changes in the value the time step at t = 0.05 and t = 0.075, i.e., at the time moments corresponding to discontinuities in the right-hand side and the coefficient of the equation. In accordance with the rule (13), the time step increases at the initial time stage. Let us decompose the correcting coefficient γ n+1 into three terms: γ f n+1 − f n , s 3 = A n+1 (y n − y n−1 ) .

Figure 2 demonstrates their influence.
The influence of the reducing error level δ on the convergence of the approximate solution is shown in Fig. 3. The approximate solution at the point x * = 0.5 is depicted in this figure. For comparison, Figure4 presents similar data that were obtained using the uniform grids in time.   Special attention should be given to the influence of the initial conditions. A typical situation is the presence of a boundary layer and this requires to use small steps at the initial time stage. For example, the behavior of the time step for our model problem with initial conditions u 0 (x) = 2x, 0 < x ≤ 0.5, 1 − 2(x − 0.5), 0.5 < x < 1, is shown in Fig. 5. Compared with Fig. 1 (smooth initial conditions), the initial time stage is calculated with essentially smaller time steps and the total number of steps is increased by more than a factor of 2. In the region outside the neighbourhood of discontinuities of the coefficients and the right-hand side, the time step is controlled first of all by the term s 3 (see Fig. 6). A more difficult situation for the numerical solution is connected with inconsistent initial and boundary conditions. Let we have u 0 (x) = 1, 0 < x < 1.
The selection of the time step for this case is shown in Fig. 7. Up to t = 0.000856 the calculation is carried out with the minimum time step τ 0 = 1 · 10 −6 . That is why the total number of time steps is 2183.