Fast explicit time integration schemes for parabolic problems in mechanics

We present a family of fast explicit time integration schemes of first, second and third order accuracy for parabolic problems in mechanics solved via standard numerical methods that have considerable higher computational efficiency versus existing explicit methods of the same order. The derivation of the new explicit schemes is inspired on the finite increment calculus (FIC) procedure used for obtaining stabilized numerical schemes in fluid and solid mechanics. The new (so‐called) explicit FIC‐Time (EFT) schemes allow considerable larger time steps than the standard first order forward Euler (FE) scheme and the second and third order Adams–Bashforth schemes. The comparison with Runge–Kutta schemes also favors the FIC‐Time schemes in terms of the limit time step size (for second order schemes) and the total number of matrix‐vector multiplications per time step (for second and third order schemes). The new first order explicit schemes have a faster convergence to steady‐state than the FE scheme. The accuracy and efficiency of the new EFT schemes are verified in examples of application to the transient heat conduction equation using the finite element method.


INTRODUCTION
Parabolic problems are common in mechanics.Examples of problems involving one single scalar variable are heat conduction, seepage in porous media, gas diffusion, magnetostatics and Reynolds film lubrication, among others. 1,2Typical multivariable parabolic problems are found in solid and fluid mechanics.][5] The explicit time integration of parabolic problems is restricted by the critical time step which governs the stability of the numerical solution (see for instance Reference 1).The value of the critical time step is proportional to the inverse of the maximum eigenvalue of the discretized system.ÑATE et al.In Reference 6 the authors presented a first order explicit time integration scheme for parabolic problems that allows larger time steps than the standard forward Euler (FE) scheme.The derivation of the new explicit scheme was based on the finite element calculus (FIC) procedure for modeling the transient heat conduction problem. 6The stability and time step size of the new explicit FIC-based scheme (hereafter termed explicit FIC-Time scheme, or EFT scheme in short) is governed by a time stabilization parameter.The examples of application in one (1D) and two (2D) dimensions presented in Reference 6 showed that considerable gains in the time step size versus the standard FE scheme were obtained with the EFT scheme, while retaining the accuracy of the transient solution.
In this work, the approach presented in Reference 6 is generalized for deriving first, second and third order EFT integration schemes for the semidiscrete parabolic equation.The new EFT schemes can be considered as a particular class of explicit linear multistep methods (LMS) 1,7 that allow time steps considerably larger than those for the standard FE and Adams-Bashforth (AB) schemes of the same order.The new EFT schemes also show advantages versus Runge-Kutta (RK) schemes in terms of the limit time step size (for second order schemes) and the number of matrix-vector multiplications per time step (for second and third order schemes).The accuracy and efficiency of the new EFT schemes are verified in three examples of application to the transient heat conduction equation using the finite element method.

THE NUMERICAL SOLUTION OF THE SEMIDISCRETE PARABOLIC EQUATION IN MECHANICS: EXPLICIT FORWARD-EULER SCHEME
The semidiscrete form of the parabolic equation in mechanics using either FEM, FVM, or FDM has the following standard form 1,2 where a(t) is the vector of unknowns of the discretized problem in space, ȧ = d dt a, f is the vector emanating from external actions to the system, and C and K are square matrices involving physical and geometrical parameters of the problem, Ω is the space domain where the parabolic equation is solved and t is the time.
Without loss of generality, in this work we will particularize the physical problem expressed by Equation (1) to transient heat conduction problems.In this context, t and Ω are respectively the time and the space domain where Equation (1) holds, a is the vector of the nodal temperatures, ( ⋅) denotes the time derivative, C is the heat capacity matrix with terms as a function of c where  is the density and c is the specific heat, K is the conductivity matrix and f is the flux vector containing contributions from the internal and external heat sources.Equation ( 1) is completed with the appropriate boundary conditions specifying the values of the temperature and the normal heat flux at the boundary, and with the initial temperature field on Ω.The expressions of the matrices and vectors of Equation (1) and the boundary conditions can be found in many text books. 1,4n this work we will assume a diagonal form of matrix C and material properties independent of time.
The explicit solution of Equation ( 1) can be written, using the FE scheme, 1,4,8 as where (⋅) n denotes values at time t = t n .From Equation (2) we can find Equation (3) shows that for a diagonal form of C the time evolution of the nodal temperatures can be found without solving a system of equations.
The stability of the FE scheme is found using standard procedures. 1,4,8A summary is given next.The solution of the discretized transient problem expressed by Equation (1) will be assumed of the form where N is the number of nodal temperatures in the discretized system and y i (t) and  i are, respectively the ith modal amplitude and the ith modal vector of the eigenvalue system In Equation ( 5)  1 , … ,  N are the eigenvalues of the discretized system in space ordered from the smallest one,  1 , to the largest one,  N .In the computation of the eigenvalues for the examples presented in the article we have disregarded the equations for the prescribed temperatures in Equation (5).
Substituting Equation ( 4) into (3) for f n = 0, and using the orthogonality properties of the modal vectors  i versus the matrices C and K, 1 we can obtain, after little algebra, the expression for the time evolution of the ith modal amplitude as and Assuming a stable solution of the type y n+1 i =  i y n i , with  i being the amplification factor for the ith mode, we can obtain from Equation (6), The stability condition, ensuring From Equations ( 7) and ( 10) we deduce the stability and time step conditions for i = N, as where Δt l is the limit time step that guarantees the stability of the explicit FE scheme.
In order to avoid non-physical oscillatory solutions, the following additional condition must be satisfied by the ith amplification factor where m is the number of low order modes that govern the solution of the transient problem.We have found that, in practice, it suffices that condition (12) is satisfied for the first mode (which typically dominates the transient solution) only, that is, 0 <  1 < 1.In the examples solved in this work this has ensured a good approximation to the correct solution in all cases.

FIC FORM OF THE SEMIDISCRETE PARABOLIC EQUATION IN MECHANICS
The differential form of the parabolic equation in mechanics is found by establishing the balance of fluxes (heat, momentum, forces, etc.) in an infinitesimal space-time domain belonging to [Ω, t] (see e.g., References 1 and 2).From the governing differential equation and the boundary conditions, the discretized form of Equation ( 1) can be found using the FEM, for instance References 1 and 4.
A higher order form of the parabolic equation can be found using a FIC procedure, as described in Reference 6.The simplest, hereafter called "FIC-Time," form of the semidiscrete parabolic equation can be written as r −  ̇r = 0, (13)   where is a residual vector of the nodal heat fluxes that vanishes for the exact solution and  is a characteristic time parameter.
For heat conduction problems Equation (3) can be found by establishing the balance of heat fluxes in a space-time domain of dimensions [dΩ, ] where dΩ = dxdy in 2D and  is a finite time interval defined as a proportion of the time step size Δt of the time integration scheme chosen.In practice,  is expressed as  = Δt, where  is a non dimensional "time stabilization parameter" that governs the stability of the time integration scheme.Clearly for  = 0 then Equation ( 13) simplifies to r = 0 which is the standard infinitesimal expression of the semidiscrete heat conduction equation, given by Equation (1).Details of the derivation of Equation ( 13) can be found in References 6,9,10.
Equation ( 13) can also be interpreted on purely mathematical grounds as a Taylor series expansion in time of the residual vector r.Appendix A shows an example of this interpretation.
][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] A general FIC-Time form of the semidiscrete parabolic equation can be written as follows where  0i ,  1j ,  2k and so forth are time stabilization parameters, r i , ̇rj , rk and so forth denote values of the residual r and its time derivatives associated with the parameters  0i ,  1j ,  2k and so forth, respectively.The sign and the terms within the residuals terms in Equation ( 15) can be formulated in a number of ways, each one leading to a different FIC-Time integration scheme, as shown in the article.In (15) the sum criterion for repeated indices applies.Again, when the time stabilization parameter are all zero, the infinitesimal form of the parabolic equation given by Equation (1) is obtained.
The time stabilization parameters  pr can be expressed in terms of the characteristic time parameters and the time increment used for solving Equation (15) as For the explicit integration schemes presented in this work the time stabilization parameters  pr are identified with a single index, or simply as .The derivation of Equations ( 15) and ( 16) is shown in Appendix A.
Equation (15) is the basis for formulating new explicit and implicit LMS time integration schemes 1,7 for solving the semidiscrete parabolic equation.In this work, Equation (15) will be the starting point for deriving first, second and third order EFT integration schemes for the semidiscrete parabolic equation that allow larger time steps than standard explicit time integration procedures.

FIRST ORDER EXPLICIT FIC-TIME INTEGRATION SCHEMES WITH LARGE TIME STEPS
A family of EFT integration schemes of first order accuracy for the semidiscrete parabolic equation can be derived starting from Equation (15).The algorithms are termed EFT1k, where the numbers 1 and k denote first order and the particular scheme, respectively.

First order EFT11 scheme
The first order EFT11 scheme was derived in Reference 6 starting from the following FIC-Time form Note that Equation ( 17) is a particular case of Equation (15).Details of the computation of the residual terms in Equation (17) are given in Reference 6.The general expression of the EFT11 scheme can be written in terms of a single stabilization parameters as 6 1 Δt

First order EFT12, EFT13, EFT14, EFT15, and EFT16 schemes
The first order EFT12-EFT16 schemes are derived starting from the following simple FIC-Time form The terms within the residuals r and r 1 can be formulated in a number of ways, leading to each of the EFT12-EFT16 schemes presented below.
EFT12 scheme EFT13 scheme EFT14 scheme EFT15 scheme EFT16 scheme The derivation of the EFT12 scheme is presented in Appendix B, as an example of the general procedure for obtaining this family of first order EFT schemes.

TA B L E 1
Relevant parameters for the six first order explicit FIC-Time (EFT) integration schemes of Sections 4.1 and 4.2.

Stability analysis and limit time step
The general expression of the different first order EFT schemes can be written in compact form as where a 1 , b 1 , c 1 , and d 1 are parameters that characterize each integration scheme.The expression of these parameters is given in Table 1.
In general, a 1 , b 1 , c 1 , and d 1 are a linear function of the time stabilization parameter .
The explicit computation of a n+1 is deduced from Equation (25) as Substituting Equation ( 4) into (25) and using the orthogonality properties of the eigenvectors  i versus the matrices C and K, we can obtain the decoupled equation for the ith modal amplitude as (for f n = 0) where p i is given by Equation (7).
The characteristic equation for the modal amplification factor,  i , is obtained as where 3]26 In either way, the three stability conditions are the following: This gives Equation ( 33) gives the condition for p i as Condition ( 34) must be satisfied for all the modes.Recalling that p i = Δt i and particularizing Equation ( 34) for the highest mode (i = N) yields the expression of the allowable time increment, as where Δt l is the limit time step ensuring a stable transient solution.
Comparing the expression of Δt l of Equation (35) with that of Equation ( 11) gives the speed-up (SUP) of the new EFT schemes versus the FE scheme (SUP FE ) as and, hence,

Critical values of the stabilization parameter and of the time step size
Following what was mentioned in the last paragraph of Section 2, in order to preserve the accuracy of the transient solution, the time step size should be chosen so that the following condition for the amplification factor of the lowest mode is satisfied Satisfaction of Equation (38) guarantees that the transient solution is free from spurious oscillations.The value of  1 can be obtained by solving the quadratic characteristic equation for  1 (Equation 28).This gives The absence of oscillatory solutions is ensured by a non-negative value of the discriminant of Equation (39), that is, Using Equation ( 29) (for i = 1), we can rewrite Equation (40) as Substituting the expressions of a 1 , b 1 , c 1 , and d 1 of each EFT scheme into Equation (41) we can obtain (for the equality case) the critical value of the stabilization parameter (denoted as  c ) that ensures no spurious oscillations along the transient solution.
On the other hand, the limit values of the stabilization parameter ( l ) must ensure that ∞ > Δt > 0. Introducing these constraints into Equation (35), we can easily deduce the following conditions for a 1 , b 1 , c 1 , and d 1 From Equations (42) the bounds of the limit stabilization parameter,  l , can be found.
Table 1 shows the value of SUP FE , the limit ( l ) and critical ( c ) values of the stabilization parameter for each of the six first order EFT11-EFT16 schemes.The table also shows the critical speed-up SUP FE c computed by substituting the expression of the critical stabilization parameter,  c , in that of SUP FE .
Table 1 shows that the value of SUP FE c is identical for the six first order EFT integration schemes here presented, which is a remarkable result.
Table 2 shows the values of the critical stabilization parameter,  c , and of SUP FE c for different eigenvalue ratios It is noticeable that SUP FE c increases as r 1 decreases (see also Equation ( 43)).Note also that time steps two orders of magnitude larger than those for the FE scheme can be achieved when r 1 ranges from 10 −4 to 10 −6 , which is a common ratio in practical problems.To be on the safe side, we recommend choosing a value for the stabilization parameter slightly smaller than the critical one.In the examples presented in the article, we have chosen  =  c , where  ranges from 0.85 to 0.985 for values of r 1 = 10 −4 − 10 −7 .This has led in all cases to a smooth transient solution computed with time steps much larger than for the FE scheme.

Remark 1. The value of SUP FE
c for the six first order EFT schemes of Sections 4.1 and 4.2 and Table 1 can be approximated as Equation ( 43) shows that the speed-up versus the FE scheme increases for larger values of the eigenvalue ratio . This will be the case for finer meshes, or for big size problems which is an outstanding result.Table 2 shows an example of this result.

TA B L E 2
Critical values of the time stabilization parameter  c and of the speed-up versus the forward Euler scheme (SUP FE c ) for different value of the eigenvalue ratio r 1 for the six first order EFT schemes of Table 1 and Sections  Remark 2. From the value of  c we can obtain the critical time step for each first order EFT scheme as Remark 3. The application of the first order EFT schemes requires knowledge of the solution at two consecutive initial times (i.e.,  0 and  1 ).The value of  1 can be obtained in a number of ways.In the examples presented in the article using the EFT12 scheme, we have computed  1 by using the FE scheme with a time step Δt = 1 10 Δt FE , where Δt FE is computed by Equation (11).Note that this implies running K FE solutions, with K = 10Δt EFT Δt FE , where Δt EFT is the time step of the first order EFT scheme chosen.Once  1 has been obtained, the transient solution carries on with the EFT scheme and a time step Δt = Δt EFT .Remark 4. The implementation complexity and the computational cost of the six first order EFT schemes of Table 1 is very similar.

Error analysis
The (absolute) error of the FE scheme can be estimated as the error for the first modal component that typically dominates the global solution.An approximate expression of the error at time t = nΔt is In Equation ( 44), y 0 is the initial solution for the first modal component, n is the number of time steps at which the error is evaluated and 44) is a good approximation of the error for small values of r 1 (say r 1 < 10 −3 ).The expression of the error for the first order EFT schemes is more cumbersome.An approximation of the error at time t is the following where Expression (45) holds for the different first order EFT schemes of Table 1.
The evaluation of Equation (45) involves choosing the time stabilization parameter  for each first order EFT scheme.This affects the time step size (see Table 1) and, consequently, the number of steps n to reach the time t for the EFT scheme chosen.
Equations ( 44) and (45) yield basically the same error for  = 0 and large values of n.For small values of n there is a small discrepancy due to the fact that in the EFT schemes we assume a quasi-exact value of  1 (see Remark 3), which is not the case in the FE scheme.
Remark 5. Negative values of the stability parameter  for the EFT11-EFT15 schemes, or smaller than 1∕4 for the EFT16 scheme yield values of SUP FE < 1, that is, the allowable time steps are smaller than for the FE scheme.Application of Equation (45) shows that the error in the numerical solution for the first order EFT schemes is smaller than that for the FE scheme in those cases.This is another remarkable property of the EFT schemes.Remark 6.For  > 0 giving values of SUP FE > 1, the error in the numerical solution of all the EFT schemes is larger than that for the FE scheme.The trade-off between speed-up in the explicit solution and the allowable numerical error should be considered for each problem.

4.6
Fast convergence to the steady-state solution An interesting application of explicit time integration schemes is their ability to obtain the steady-state numerical solution of a mechanical problem (i.e., that of the system Ka = f) without solving a system of equations.][29][30][31][32] Each time increment can be considered a step towards the steady-state solution, 32 which is the target result, irrespective of the results along the transient solution.This, of course, only typically holds for linear (or linearized) problems where K and f are independent of a.
It can be shown that the equation governing the error between the steady-state solution and the transient solution at each time step for each mode has an identical form to Equation (28) for the modal amplification factor, simply substituting  i by  i , where  i is the error for the ith mode.
For instance, the modal error versus the steady-state solution per time step for the FE scheme is obtained as In practice, following what has been mentioned earlier, the error of the overall transient solution is governed by that of the first mode.Particularizing Equation (48) for i = 1 we can obtain an approximate expression of the error towards the steady-state solution after M time steps for the FE scheme as Clearly as r 1 is typically very small for practical problems, the convergence of the FE scheme towards steady-state is typically very slow.
All the new first order EFT schemes of Table 1 converge considerably faster to steady-state than the FE scheme.It can be shown that the faster convergence is obtained for the critical value of the stabilization parameter,  c .The error per time step in this case is obtained for the six EFT schemes of Table 1 as The derivation of Equation (50) for the EFT12 scheme is given in Appendix B.
Table 1 shows the speed-up in the convergence toward the steady-state solution versus the FE scheme (SUP FE c ).We note again that this value is the same for the six first order EFT schemes here presented.
We highlight that expression (50) is identical to the error towards steady-state in explicit schemes for second order dynamic problems with optimal Rayleigh damping. 29,31This opens the door for deriving enhanced time relaxation schemes in computational mechanics via the FIC-Time procedure.
Remark 7. If the steady-solution is sought, the expression of the diagonal matrix C can be chosen in a convenient arbitrary form so as to maximize the value of r 1 (Equation 49).A usual choice in time relaxation schemes is C ii = K ii .

Second order explicit EFT21, EFT22, EFT23, and EFT24 schemes
The general expression of the second order EFT schemes is obtained starting from the following FIC-Time form Equation ( 51) is a particular case of Equation (15).The terms within the residuals r, r 1 , and ̇r1 are formulated in a number of ways leading to different second order EFT schemes.We present below four members of this family, termed EFT21, EFT22, EFT23, and EFT24.
EFT21 scheme ] EFT22 scheme EFT23 scheme EFT24 scheme Grouping terms, the form of Equation ( 25) can be found for each of the above second order EFT schemes with the parameters a 1 , b 1 , c 1 , and d 1 given in Table 3.

Stability analysis and limit time step
The general form of the four second order explicit EFT schemes given in the previous section can be written as The explicit computation of a n+1 is deduced from Equation (56) as The characteristic equation for the modal amplification factor for the ith mode,  i , is obtained as TA B L E 3 Relevant parameters for the four second order EFT schemes of Section 5.1.
Note: Speed-up (SUP) versus the second order Adams-Bashforth (AB2) and Runge-Kutta (RK2) schemes. where The four stability conditions for the coefficients A, B, C, D ensuring −1 <  i < 1 are in this case (a) The conditions (60) can be derived via the Routh-Hurwitz criterion. 1,3,26The derivation is more cumbersome in this case as the characteristic equation, Equation (58), is cubic.Note that condition (60d) is nonlinear, which complicates the derivation of stable second order EFT schemes.
The limit time step is typically governed by condition (60c).Substituting A, B, C, and D from (59) into (60c) gives From (61) we deduce The limit time step is obtained recalling the expression of p i (Equation 7) and making i = N, that is, Table 3 shows the value of a 2 , b 2 , c 2 , d 2 , e 2 , and f 2 for the four second order EFT schemes given in Section 5.1.
Remark 8.The limit time step size of the second order EFT22, EFT23, and EFT24 schemes is 2  N (i.e., as for the FE scheme, but with higher accuracy).For the EFT21 scheme, a considerably larger limit step can be used as  1 approaches the unity.This makes the second order EFT21 scheme the better candidate for solving parabolic problems using larger time steps than for standard second order explicit schemes.The time step size is, however, smaller than for the first order EFT schemes presented in Section 4.
Remark 9.If f 2 = 0 in Equation ( 59), then the four stability condition of Equation (60) reduce to the three conditions expressed by Equations ( 30), (32), and (33) for first order schemes, which simplifies the stability analysis of the EFT22, EFT23, and EFT24 schemes.
Remark 10.Application of the second order EFT schemes requires knowledge of the initial solution at three consecutive times (i.e.,  0 ,  1 , and  2 ).Computation of  1 and  2 in this work has been performed in a similar way as explained in Remark 3 for the first order EFT schemes.

THIRD ORDER EXPLICIT FIC-TIME INTEGRATION SCHEME
The matrix expression of the third order EFT integration scheme here presented, termed EFT31, is written as Third order accuracy is obtained by imposing the following condition on the stabilization parameter  2 Substituting Equation (65) into Equation (64) and grouping terms leads to the following compact expression for the third order EFT31 scheme ) . (67) The characteristic equation is cubic and the stability analysis follows the same lines expressed by Equations ( 58)-(60d).
The limit time step size can be found as with the following restrictions for  1 and  2 , It can be found that The third order EFT31 scheme has the same limit time step as the FE scheme, but it has higher accuracy.

COMPARISON OF THE SECOND AND THIRD ORDER EXPLICIT FIC-TIME SCHEMES WITH THE EXPLICIT ADAMS-BASHFORTH AND RUNGE-KUTTA SCHEMES
For comparison purposes we will consider the second and third order explicit AB and RK schemes.For simplicity we will assume f = 0 in all cases.

Second and third order Adams-Bashforth schemes (AB2 and AB3)
The explicit AB schemes are a particular class of explicit LMS methods. 1,7We present below the AB2 and AB3 schemes.

AB2 scheme
It is expressed as The limit time step is given by AB3 scheme The limit time step is

7.2
Second and third order Runge-Kutta schemes (RK2 and RK3) RK schemes are multi-stage methods. 33,34The explicit nth RK scheme involves n multiplications of the system matrix with the vector of nodal unknowns at each time step.
RK2 scheme with The limit time step is RK3 scheme with ) , The limit time step is In Equations ( 76) and (79) K is given by Equation (71b).

Comparison of the second and third order EFT schemes with AB and RK schemes
Table 3 shows the speed-ups of the second order EFT schemes versus the AB2 and RK2 schemes (termed SUP AB2 and SUP RK2 , respectively) in terms of the stabilization parameters for each case.We note the following: -The EFT21 scheme allows considerably larger time steps than the AB2 and RK2 schemes.
-The EFT22, EFT23, and EFT24 schemes allow time steps about twice larger than the AB2 scheme and about the same value than the RK2 scheme.
The third order EFT31 scheme has the following speed-up versus the third-order Adams-Bashforth (AB3) and Runge-Kutta (RK3) schemes Expressions (81) show that the third order EFT31 scheme allows time steps some four times larger than the AB3 scheme for the same number of matrix-vector multiplications per time step.On the other hand, the EFT31 scheme is slightly less competitive than the RK3 scheme in terms of limit step size, although it requires one third of the number of matrix-vector multiplications per time step, as explained in the next section.

COMPACT FORM OF THE EFT AND AB SCHEMES AND NUMBER OF OPERATIONS VERSUS RK SCHEMES
As it is typical in explicit LMS methods, 7 the different EFT and AB integration schemes can be expressed in a compact form as (for f = 0) with K given by Equation (71a), and where  0 ,  1  2 are parameters of each explicit LMS scheme.Equation ( 82) clearly shows that both the EFT and AB methods require storing the "historical data pool."However, they only require a single multiplication of a matrix, K, by a vector, ân , at each time increment for any order of the scheme chosen.
On the other hand, RK schemes, as a class of multi-stage methods, require a number of matrix-vector multiplications equal to the order of the RK scheme, that is, two matrix-vector multiplications for RK2, three for RK3 and so forth.
This fact makes the EFT schemes competitive in term of the number of operations versus the RK schemes of the same order.We have verified this property for selected second and third order EFT schemes in the example shown in Section 10.3.

SUPER-TIME STEPPING ACCELERATION OF THE FORWARD EULER SCHEME
The super-time stepping (STS) scheme for accelerating the FE scheme was introduced, as early as in 1978, by Gentzsch and Schlüter. 35More details were given in Reference 36.The method was applied in Reference 37 within the context of multigrid procedures for elliptic problems.The STS scheme was brought back to the attention of the computational community in 1996 by Alexiades et al. 38 The STS scheme basically relaxes the requirement of stability at the end of each time step to requiring it at the end of a cycle of N s sub-steps.
The STS scheme can be considered as a class of Runge-Kutta-Chebyshev methods (see Reference 38 and references therein).It can be formulated for the solution of Equation ( 25) as where j refers to the jth of the N s supersteps dividing the time step Δt chosen, that is, Δt = ∑ N j=1  j .The limit value of  j is given

𝜆
) as In finite element practice,  is estimated by taking  N = max  e where max  e denotes the largest eigenvalue of all the individual elements in a finite element mesh, 1,3,4 and  1 is computed, using a coarse mesh for a given problem.The STS scheme is first order accurate.
The STS scheme allows limit values of Δt considerably larger than for the FE scheme.The prize to be paid is the need to compute N s intermediate solutions, each one requiring a matrix-vector operation, as in RK schemes.The actual SUP FE value of the STS scheme should be computed by dividing the number of solutions steps by N s .In some problems, however, the STS scheme can beat the FE scheme in computational efficiency.
It should be finally noted that the intermediate solutions, a n+j have typically no physical meaning in the STS scheme, which, somehow, is a drawback of this method.

Transient heat conduction in square domain with internal circle
The example chosen to show the performance of selected first and second order EFT schemes is the time evolution of the temperature in a square prismatic domain of unit thickness with an internal circle (Figure 1).We have solved the transient problem with the following first order explicit schemes: FE, EFT12 and STS.The problem has also been solved with the second order EFT21 scheme.
The material properties are the following, The distributed heat source Q is zero.The initial temperature at the internal points of the domain is prescribed to zero, with exception of the values at the mid-length of the lateral sides where the initial temperature is prescribed (and kept fixed) at 100 and 0 • C, as shown in Figure 1.A zero normal flux is specified at all the other boundaries.
The problem has been analyzed with the four unstructured meshes of of 870, 3432, 13,732, and 59,924 3-noded linear triangles shown in Figure 2. The number of nodes, the minimum and maximum eigenvalues and the value of the critical stabilization parameter for the EFT12 scheme for each mesh are given in the figure .The curves in Figure 3 display the evolution of the temperature from the initial value to steady-state at the points A and B shown in Figure 1.We discuss below the results obtained with each explicit scheme.

F I G U R E 1
Transient heat conduction in square domain with internal hole.Geometry, material properties, boundary and initial temperature conditions at all the internal points of the domain.Temperature contours at steady-state obtained with the forward Euler scheme and a mesh of 3432 3-noded triangles.

F I G U R E 2
Four unstructured meshes of 870, 3,432, 13,732, and 59,924 3-noded linear triangles used for solving the problem of Figure 1.
Values of the larger and smaller eigenvalues for each mesh and of the critical stabilization parameter for the first order EFT12 scheme.

Forward Euler (FE) scheme
Figure 3 shows the time step used for the FE scheme for each mesh.The time step for each case has been computed with Equation (11).The curves for the temperature evolution at points A and B obtained with the FE scheme for each mesh (Figure 3) are taken as the reference for comparison with the other explicit schemes chosen, as discussed below.
Figure 1 shows the steady-state distribution of the temperature over the analysis domain obtained with mesh 2 of 3432 elements and the FE scheme.The steady-state solution is very similar for the four meshes considered.

First order EFT12 scheme
Figure 3 shows the value of the stabilization parameter  used for each mesh in the computations with the EFT12 scheme.We have verified that the values of  are slightly smaller than the critical values shown in Figure 2. The values chosen for  yield accurate transient solutions in all points of the domain.The actual time steps used for each analysis are shown in Figure 3.
Excellent agreement between the EFT12 and the FE results for the temperature evolution at points A and B of Figure 1 is obtained in all cases.Similar good agreement is found in the rest of the analysis domain, as well as for the steady state solution, that basically coincides with that shown in Figure 1.
F I G U R E 3 Time evolution of the temperature at points A and B of Figure 1 using the forward Euler (FE), EFT12, EFT21, and STS schemes for the four meshes of Figure 1.The solution for all the schemes is practically coincident.
The values of SUP FE for the ETF12 scheme are shown for each mesh.These values show that the gains in time step size versus the FE scheme range from 11.3 times for the coarser mesh, to 83.46 times for the finer one.The advantages of the EFT12 scheme versus the FE scheme in terms of time step size are clearly evidenced.
Figure 4 shows the temperature evolution at points A and B obtained with the EFT12 scheme for the mesh of 3,432 elements when the time step is smaller (continuous line) and greater (dotted line) than the critical value  c .The correct steady state solution is obtained in both cases, as expected, despite the fact that the transient solution has no physical meaning for time steps larger than the critical one.This shows the importance of using time steps equal or smaller than the critical value for the first order EFT schemes, if meaningful transient solutions are sought.
The above conclusions apply to the rest of the first order EFT schemes of Table 1.

Second order EFT21 scheme
Figure 3 shows the results for the temperature evolution at points A and B obtained with the second order EFT21 scheme.Excellent agreement between the solutions for the FE and EFT21 schemes is found for all meshes.We note that gains in time step size versus the FE scheme are smaller than for the EFT12 scheme (3.97 times for the coarser mesh to 11 times for the finer one).However, the time steps of the EFT21 scheme for all cases are still considerably larger than for the AB2 scheme (see Equation ( 72)) with equivalent good results in both cases.In any case, the efficiency of EFT21 scheme shows clearlier in terms of its accuracy versus the number of algebraic operations per time step, in addition to the time step size (see Section 10.3).

STS scheme
Figure 3 displays the results obtained with the STS scheme, for the sake of comparison.Larger time steps than for the FE scheme are obtained.Note, however, that the actual time step size should be divided by the number of sub-steps (N s ), each one implying a matrix-vector multiplication.The SUP FE value shown in Figure 3 is the result of this quotient.In the whole, the EFT12 and EFT21 schemes beat the STS scheme in terms of computational efficiency in this example for all the meshes considered.

Convergence towards steady-state
We have evaluated the efficiency of the explicit time integration schemes here considered for reaching the steady-state solution of the heat conduction problem of Figure and  = 10 −6 ).For comparison purposes, the steady-state solution was also obtained by solving the system of equations Ka = f directly.
Figure 5 shows the error versus the steady-state solution at point A obtained using the four explicit integration schemes chosen: EF, ETF12 (with  =  c ), ETF21 (with the values of  given in Figure 3) and STS for the four meshes considered.The faster convergence in all cases is obtained with the EFT12 scheme using the critical stabilization parameter  c .The speed-up in the convergence toward the steady-state solution for the EFT12 scheme versus the FE scheme is equal to SUP FE c (see Section 4.6 and Remark 1).(with  =  c ), EFT21 (with  1 given in Figure 3) and STS schemes.The faster convergence in all cases is obtained with the EFT12 scheme using the critical stabilization parameter  c .
3, and 4 of Figure 2, respectively.Note that, as expected, the value of SUP FE c increases for finer meshes for which the ratio increases (Remark 1).

Discontinuous sinusoidal distributed heat source pulse acting within a rectangular domain
The geometry of the problem, the material properties and the boundary and initial conditions are shown in Figure 6.The temperature is prescribed to • C at the left vertical side.A zero normal flux is assumed to the other three sides.The discontinuous sinusoidal heat source pulse acting in the internal circular domain (marked in blue) is shown in Figure 6.
The period of the heat source pulse is 10 s.The heat pulse features constraint the maximum time step that can be applied in implicit schemes, for the sake of accuracy.
The numerical study has been performed using the following integration schemes: -First order explicit schemes: FE, EFT12 and STS.
-Backward Euler (BE) implicit scheme.The system of algebraic equations for each implicit time step has been solved iteratively with the standard conjugate gradient (CG) method and directly with the Cholesky method.The "exact" solution of this problem was taken as the one obtained with a mesh of 1,150,755 3-noded triangles using the BE scheme with Δt = 0.01 s.
Figure 7 shows a table with the time step size for each numerical scheme, the total number of steps to reach a solution at a time of 1000 seconds and the total CPU time for each scheme for a coarse mesh of 452 3-noded triangles.The figure displays the numerical results for the time evolution of the temperature at the points A and B of Figure 6.The plotted curves shown are quite undistinguishable for all the schemes considered.Figure 7 presents the speed-up (SUP) in terms of number of solution steps and total CPU time, taking as a reference the values for the FE scheme.It is noticeable that the EFT12 scheme has a CPU cost twice smaller than the FE scheme.The figures of the SUP in CPU time also indicate that the implicit BE scheme (using the CG method) is not competitive in this case as it requires one order of magnitude larger total CPU time than the EFT12 scheme.The Cholesky method for solving the equations in the BE scheme has a CPU considerably larger than for the CG method, and it has been disregarded in subsequent solutions of this problem.
Figure 8 shows the error in the results along the transient solution at point A for the different integration schemes considered.The maximum error corresponds to the implicit BE scheme using the CG method.
A similar study is presented in Figures 9 and 10 for a fine mesh of 184,992 3-noded triangles.The EFT12 scheme is again the more competitive solution method in terms of total CPU time, with gains of 41.80 times versus the FE scheme and of about four times versus the implicit BE scheme.
The EFT12 scheme also yields an smaller error than the BE scheme in the temperature values along time at point A (Figure 10).
Figure 11 shows the SUP values in terms of gains in CPU time for the four time integration schemes considered, versus the number of degrees of freedom (DOF) for different meshes of 3-noded triangles used for solving this problem.The efficiency of the EFT12 scheme is clearly evidenced.
It is interesting that the STS scheme is quite competitive for solving this example despite that its efficiency is smaller than that for the EFT12 scheme for all meshes.

F I G U R E 8
Discontinuous sinusoidal heat source pulse on rectangular domain.Error in the time distribution of the temperature at points A and B of Figure 6 for the mesh of 452 3-noded triangles using the FE, EFT12, STS, and BE (using CG and Cholesky methods) schemes.

F I G U R E 9
Discontinuous sinusoidal heat source pulse on rectangular domain.Upper table: Time step size for the FE, EFT12, STS and backward Euler schemes; total number of steps to reach a solution at 1000 s. and CPU time for each scheme using 184,992 3-noded triangles.Intermediate curves: Numerical results for the time evolution of the temperature at the points A and B of Figure 6.Lower graphs: speed-up (SUP) in terms of number of solution steps and of the total CPU time, taking as a reference the values for the forward Euler (FE) scheme.

F I G U R E 10
Discontinuous sinusoidal heat source pulse on rectangular domain.Error in the time distribution of the temperature at points A and B of Figure 6 for the mesh of 184,992 3-noded triangles using the FE, EFT12, STS, and BE (with the CG method) schemes.The efficiency of the implicit solution can be improved by using a pre-conditioned CG method.The comparison with the new explicit scheme will be the object of future studies.

Transient heat conduction along a bar: Error analysis of the EFT12, EFT21, and EFT31 schemes
The simple problem considered for performing an error analysis of the EFT12 and EFT21 schemes is the time evolution of the temperature in a bar. Figure 12 shows the geometry of the bar, the mechanical properties and the boundary and initial conditions for the temperature along the bar.The later has been chosen so as to ensure a smooth transition of the temperature at all points of the bar from the initial condition towards the steady-state solution, characterized by a linear distribution of the temperature between the two bar ends.The problem has been solved with a uniform mesh of 40 2-noded elements.
The reference temperature T is obtained using the 4th order Adams-Bashforth scheme 7 with the same mesh and a very small step Δt = Δt E 2 15 , where Δt E denotes here the limit time step of the FE scheme (Equation 11).At a given simulation time t, we compute the root mean square error (RMSE) of the numerical solution as In Equation (87) T ref (x) and T num (x) are, respectively, the values of the reference temperature and the temperature computed numerically at the x coordinate.
Figure 13 presents the convergence study of two first order schemes: the FE scheme with the limit time step, Δt E , and the EFT12 scheme using different time steps as a proportion of Δt E .The curves show the error obtained at different times.Results for the EFT12 scheme for time steps larger than Δt E (i.e., for SUP FE > 1) have a bigger error than for the FE scheme.These errors, however, can be acceptable for practical purposes.On the other hand, the EFT12 scheme with time steps smaller than Δt E , has smaller errors than the FE scheme.Figure 14 (left) shows a similar set of curves for the second order EFT21, AB2 and RK2 schemes.The time steps are compared with that of the FE scheme.Results for the EFT21 scheme have a smaller error than for the AB2 scheme for larger time steps.On the other hand, the EFT21 scheme results have a slighter larger error than for the RK2 scheme for similar time step sizes.
The results in Figure 14 (right) are more conclusive of the competitiveness of the EFT21 scheme versus the AB2 and RK2 schemes.The EFT21 scheme requires less number of operations to reach a given error versus the AB2 and RK2 schemes in all the cases here considered.
Figure 15 (left and right) present a similar comparison of the error obtained with the third order EFT31, AB3 and RK3 schemes.The same conclusions on the good performance of the EFT31 scheme in terms the number of operations for obtaining a certain error, as mentioned in the previous paragraph for the EFT21 scheme, are found.

CONCLUSIONS
A family of new explicit LMS time integration schemes for the semidiscrete parabolic equation has been derived starting from the general FIC form of this equation.The new explicit integration methods have been termed EFT schemes.The performance of the six first order EFT integration schemes here studied (EFT11-EFT16) is quite identical.All the six first order EFT schemes allow much larger time steps than the FE scheme.
We have also presented four members of the family of second order EFT schemes (EFT21-EFT24) and a member of the third order family of this class of new explicit schemes (EFT31).The second order EFT21 scheme and the third order EFT31 scheme have been studied in some detail.The conclusions are that both EFT schemes allow larger time steps than the corresponding explicit second and third order AB and RK schemes.Also both the EFT21 and EFT31 schemes show a good behavior in terms of the total number of operations to obtain a prescribed error in the transient numerical solution versus the AB and RK schemes of the same order.
New time integration schemes (both explicit and implicit) can be derived starting from the general FIC-Time form of the semidiscrete parabolic equation presented in this work.The competitiveness of these new integration schemes versus existing one has to be assessed in detail. with The critical value of  1 is obtained from Equation (B6) for  =  c as The critical speed-up versus the FE scheme is obtained as The expression of the error towards the steady-state solution per time step, , coincides with Equation (B5a), substituting  1 by  1 .The minimum error is given by an expression identical to Equation (B9), that is, Substituting  c from of (B8a) into (B11) gives From the expressions of R 1 and G 1 of Equations ( B5a) and (B8b) we can deduce For small values of r 1 , then Substituting (B14) into (B12) gives Equation (B15) coincides with the error towards the steady-state solution in explicit time relaxation schemes, based on solving the second order dynamic equation with a central difference scheme and optimal Rayleigh damping. 31

F I G U R E 4
Heat conduction in square domain with internal hole.Time evolution of the temperature at points A and B obtained the mesh of 3,432 elements and the EFT12 scheme when the time step is smaller (continuous line) and greater (dotted line) than the critical value  c shown in Figure2.

F I G U R E 5
Heat conduction in square domain with internal hole.Error versus the steady-state solution at point A for the FE, ETF12

F I G U R E 6
Discontinuous sinusoidal heat source pulse acting on a circular region within an internal rectangular domain, marked in blue.Geometry, material properties and boundary and initial conditions.Time distribution of the heat pulse.

F I G U R E 7
Discontinuous sinusoidal heat source pulse on rectangular domain.Upper box: Time step size for the FE, EFT12, STS and backward Euler (BE) schemes; total number of steps to reach a solution at 1000 s. and CPU time for each scheme using 452 3-noded triangles.Intermediate curves: Numerical results for the time evolution of the temperature at points A and B of Figure6.Lower graphs: speed-up (SUP) in terms of number of solution steps and total CPU time, taking as a reference the values for the forward Euler (FE) scheme.

F
I G U R E 11 Discontinuous sinusoidal heat source pulse on rectangular domain.SUP values in terms of CPU time for the FE, EFT12, STS, and BE schemes used for solving the problem with different meshes of 3-noded triangles characterized by the total number of degrees of freedom (DOF).

F I G U R E 12
Time evolution of the temperature in a bar.Geometry, mechanical properties and boundary and initial conditions for the temperature (in Kelvin) along the bar.Δt E is the limit time step of the forward Euler scheme.F I G U R E 13 Time evolution of the temperature in a bar.Value of the RMSE along time for different time step sizes for the first order forward Euler (FE) and EFT12 schemes.

F I G U R E 14
Time evolution of the temperature in a bar.Value of the RMSE along time for different time step sizes for the second order EFT21, AB2, and RK2 schemes.Δt E denotes the limit time step of the forward Euler scheme.F I G U R E 15 Time evolution of the temperature in a bar.Value of the RMSE along time for different time step sizes for the third order EFT31, AB3, and RK3 schemes.Δt E denotes the limit time step of the forward Euler scheme.
where  ref is a reference temperature and  is a very small number (here we have taken  ref = 100 This gives speed-up values of 25.31, 45.45, 86,95, and 194,17 for the meshes 1, 2,