High-Order Semi-implicit Schemes for Evolutionary Partial Differential Equations with Higher Order Derivatives

The aim of this work is to apply a semi-implicit (SI) strategy in an implicit-explicit (IMEX) Runge–Kutta (RK) setting introduced in Boscarino et al. (J Sci Comput 68:975–1001, 2016) to a sequence of 1D time-dependent partial differential equations (PDEs) with high order spatial derivatives. This strategy gives a great flexibility to treat these equations, and allows the construction of simple linearly implicit schemes without any Newton’s iteration. Furthermore, the SI IMEX-RK schemes so designed does not need any severe time step restriction that usually one has using explicit methods for the stability, i.e. Δt=O(Δtk)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t = {\mathcal {O}}(\Delta t^k)$$\end{document} for the kth (k≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k \ge 2$$\end{document}) order PDEs. For the space discretization, this strategy is combined with finite difference schemes. We illustrate the effectiveness of the schemes with many applications to dissipative, dispersive and biharmonic-type equations. Numerical experiments show that the proposed schemes are stable and can achieve optimal orders of accuracy.


Introduction
Many time-dependent PDEs which arise in physics or engineering involve the computation of high order spatial derivatives. In this paper, we are interested in proposing a semi-implicit (SI) approach with an IMEX RK setting for solving several types of PDEs with high order spatial derivatives.
We consider a sequence of such PDEs with increasingly higher order derivatives. To simplify the presentation, the PDEs examples below are only one-dimensional equations.
• The second order diffusion problem where a(u) ≥ 0 is smooth and bounded and it is a PDE with second order derivatives. Many PDE of the form (1) which arise in physics or engineering, usually involve the computation of nonlinear diffusion terms, such as: the miscible displacement in porous media [1] which is widely used in the exploration of underground water, oil, and gas, the carburizing model [2] which is derived in the chemical heat treatment in mechanical industry, the high-field model in semiconductor device simulations [3,4], and so on.
In this paper we also consider one-dimensional version of the convection-diffusion equation • The dispersive equation [5,6] with third derivatives where f (u), r (u) and g(u) are arbitrary (smooth) functions. The Korteweg-de Vries (KdV) equation [7] which is widely studied in fluid dynamics and plasma physics, is a special case of Eq. (3) for the choice of the functions f (u) = u 2 , g(u) = u and r (u) = u. The KdV-type equations play an important role in the long-term evolution of initial data [8], are often used to model the propagation of waves in a variety of nonlinear and dispersive media [9]. Another choice of the functions f (u) = u 3 , r (u) = u 2 and g(u) = u/2, gives the so called general KdV equation [6] u t + (u 3 Equation (4) is known to have compacton solutions of the form: Furthermore, Eq. (4) is a particular case of the nonlinear dispersive equation [10] u t + (u m ) x + (u(u n x x )) x = 0, m > 1, m = n + 1, with m = 3 and n = 2. In the numerical tests section we will consider Eq. (6) with m = 3, n = 2 and m = 2, n = 1. Note that in general, the prototype of nonlinear dispersive equations is the K (m; n) Eq. (6), introduced by Rosenau and Hyman in [11]. For certain values of m and n, the K (m; n) equation has solitary waves which are compactly supported. These structures, the so-called compactons, have several things in common with soliton solutions of the Korteweg-de Vries (KdV) equation where a nonlinear dispersion term replaces the linear dispersion term in the KdV equation. • The fourth order diffusion equation is a special biharmonic-type equation, where the nonlinearity could be more general but we just present (7) as an example. In this paper we will also concentrate on the one dimensional case biharmonic type equation The fourth order diffusion problem has wide applications in the modeling of thin beams and plates, strain gradient elasticity, and phase separation in binary mixtures [12].
For all these equations suitable initial conditions and boundary conditions will be set. It is well known that the time discretization is a very important issue for time dependent PDEs. For the k-th (k ≥ 2) order PDEs, explicit methods always suffer from stringent and severe time step restriction, i.e., t = O( x k ), for the stability where t is the time step and x is the mesh size. This time step is too small, resulting in excessive computational cost and rendering the explicit schemes impractical. On the other hand, implicit methods can overcome the drastic time step size restriction imposed by the stability condition for explicit schemes. However, nonlinear algebraic system must be solved (e.g. by Newton iteration) at each time step.
To cope with both the shortcomings of the explicit and implicit methods, one possible approach is to use implicit-explicit (IMEX) methods. These methods have been proposed and studied by many authors, for example [13][14][15][16]. IMEX schemes have been successfully applied to various problems, such as convection-diffusion-reaction systems [13,16,17], hyperbolic systems with relaxation [15,18], and collisional kinetic equations [19].
In the case of evolutionary partial differential equations with higher-order derivatives, IMEX methods can be used to treat the different derivative terms differently. Specifically, the higher-order derivative terms are treated implicitly, while the rest of the terms are treated explicitly. Using IMEX methods can help to alleviate the stringent time step restriction and reduce the difficulty of solving algebraic equations when the higher-order derivative terms are linear. However, for equations with nonlinear high derivative terms, IMEX methods may be much more expensive than explicit methods because nonlinear algebraic system must be solved (e.g. by Newton iteration) at each time step. In order to overcome this constraint in [20] the authors introduced the explicit-implicit-null (EIN) time-marching method. This method consists to add and subtract a sufficiently large linear highest derivative term on one side of the considered equation. After that, the linear highest derivative term is treated implicitly, while the remaining term is treated explicitly using an IMEX R-K setting.
We point out that the technique by adding and subtracting the same quantity on one side of the considered equation has been used in several papers for treating hyperbolic system with a stiff source term. For example in the paper [21], the authors for solving numerically the Bolzmann kinetic equation, due to the nonlinear stiff collision (source) term, a penalization technique have been proposed. This technique consists to penalize the Bolzmann collision operator by adding and subtracting the BGK operator and to discretize in time by an implicit-explicit approach so that the general Boltzmann operator can be handled as easily as the much simpler BGK operator. Another example has been introduced in the context of hyperbolic systems with stiff relaxation in the so-called diffusion limit. Here, standard IMEX-RK schemes are useless because such schemes become consistent discretization of the limit diffusive relaxed system, but suffer from the parabolic CFL restriction. In [22][23][24] this drawback was overcame by a penalization method, which was based on the addition of two opposite diffusive terms in the limit of large stiffness so that the limit scheme becomes an implicit scheme, unconditionally stable, for the limit diffusive relaxed system.
Generally if the leading high order spatial derivative term in (2), (3) and (8) is linear, e.g. for the Eq. (2): Eq. (3): and Eq. (8): with d > 0, a classical IMEX-RK method can be used, where the convection term is treat explicitly and the leading high order term implicitly. From a stability point of view we know that for the convection-diffusion and convection-biharmonic type equations (9), (11) if d is large enough the IMEX scheme is unconditionally stable, i.e., there exist a positive constant τ 0 independent of the spatial mesh size x, such that if t ≤ τ 0 , then the solution is strongly L 2 stable, [17,25,26]. If d is very small and the mesh size x is not too small, then the scheme is stable under the usual CFL condition for the explicit time stepping of the convection equation, i.e. t ≤ C x, where C is the Courant number, [13,16,17,25,26]. Basically, in practice we can use the stability condition t ≤ max(C x, τ 0 ). For the convection-dispersion type equation (10), the IMEX scheme is strongly stable under the classical CFL condition t ≤ C x, for hyperbolic conservation laws, see [20,27]. Although this analysis is only performed on the linear equations (9)-(10)-(11) containing the highest derivatives, numerical experiments show that the stability criterions appear to be also valid for nonlinear equations (2), (3), and (8) as shown by the numerical experiments in Sect. 3.
In this paper we discuss an alternative approach to [20] for solving time dependent PDEs with high order spatial derivatives by using high order IMEX-RK methods in a much more general context than usually found in the literature, obtaining very effective schemes. The basic idea is to propose a new strategy, called semi-implicit (SI), based on IMEX Runge-Kutta methods (SI-IMEX-RK), for the construction of a class of schemes for the solution of Eqs.
The above strategy is really convenient and useful in the case in which we have a linearly implicit evaluation for the unknown variable in the term involving high order spatial derivatives, as in Eqs. (1), (2), (6), (7) and (8). The linearly implicit evaluation is the key to the methods working for our problems. However, in other cases for time-dependent PDEs with fully nonlinear high order spatial derivatives, as for example Eq. (3), the approach proposed in [20] is more suitable.
In order to apply the SI-IMEX-RK idea, we take the nonlinear diffusion equation (1) as an example to introduce the approach in detail. Assume that the semi-discrete formulation of (1) can be written as where is the approximate solution at spatial position x i for i = 1, . . . , M, i.e., U i (t) ≈ (u(x i , t), for i = 1, . . . , M, and uniform grid spacing x i = x i+1 −x i . B ∈ R M×M is a tridiagonal matrix arising from the discretization of (a(u)u x ) x . Here we emphasize that the matrix B inherits its discontinuous dependence on U from that of a(u) on u. The SI-IMEX-RK approach is based in choosing the time discretization by an implicit and explicit RK scheme, respectively, of an IMEX pair of schemes accordingly. Roughly speaking, in the product B(U )U the implicit treatment is applied only to that second factor U , while the nonlinear term B(U ) is treated explicitly. This approach does not require solutions for nonlinear systems since the new methods require only solving a discretized diffusion equation with a linear diffusion term in which the matrix B(U ) is given.
An advantage of this approach with respect to [20] is that we can use different types of IMEX-RK schemes already exiting in the literature. In [20] the drawback in the technique of adding and subtracting the same term on one side of the considered equation is that the constant that stabilizes the scheme depends on the specific type of IMEX RK method selected. We will show several examples of equations of the form (1), (2), (6), (7) and (8) that can be efficiently solved with the SI-IMEX-RK approach, choosing different IMEX RK tableaux.
In this paper, we coupled high order finite difference schemes with high-order SI-IMEX-RK time discretization for solving diffusion (1), (2), dispersive (6) and fourth order diffusion equations (7), (8), respectively. We choose the finite difference schemes to discretize the space because of its simplicity in design and coding. However, other type of space discretization can be considered as local discontinuous Galerkin schemes [31] with application to various high order PDEs, [5,6,26,[32][33][34][35].
The organization of the paper is as follows. In Sect. 2 we present the numerical scheme by the combination of a high order semi-implicit time discretization with a suitable spatial discretization. Section 3 shows a series of numerical experiments to evaluate the performance of the proposed SI-IMEX-RK approach. Finally concluding remarks are given.

The Numerical Scheme
In this section, the strategy of numerical discretization is different from the traditional methodof-lines approach, since we first perform time discretization, after which we apply a suitable space discretization. The feature of the scheme is the following: we design a SI-IMEX-RK time discretization strategy, so that the scheme is stable with no stringent time stepping constraint and can be implemented in a semi-implicit manner [23,28,29] to enable effective and efficient numerical implementations; meanwhile we adopt high order spatial discretization strategy, so that the scheme can successfully solve the Eqs. (2), (6), (7) and (8), respectively. Next we will present a first order semi-implicit temporal discretization applied to Eq. (12), and then discuss high order extensions afterwards.

First Order Semi-implicit Temporal Discretization
It is known that in the one-dimensional case, discretizing Eq. (1) by coupling the classical explicit first order in time and second order in space method, the numerical scheme is conditionally stable, with the size of the time step: t = x 2 / 2 max U n j (B(U n j )) . Now, we couple the simplest first order semi-implicit scheme (forward-backward Euler method) with centered scheme and show that the scheme is unconditionally stable. First we approximate in space (1) at points halfway between the grid points x j , with a j+1/2 = a(U j+1 (t)) + a(U j (t)) 2 , and the analogous approximation at x j−1/2 . Differencing these then we get a centered approximation to (1) at the grid point x j Then equation (1) can be written in the form (12), where the matrix B(U (t)) has the form ⎡ Now applying the first order SI scheme to (12) we get where the numerical solution is given by Then we state that this scheme is unconditionally stable, in the sense that ||U n || ≤ ||U 0 || for all n and for any positive time-step t > 0, which is guaranteed provided that the spectral radius In order to guarantee inequality (15) we have to prove that the matrix B(U n ) is nonsingular and negative definite. The matrix (13) has the advantage of being symmetric, diagonally dominant. Moreover, since a j+1/2 , a j−1/2 > 0 for j = 1, . . . , M (as a > 0), the matrix (13) is irreducible 1 and then it can be shown that the matrix −B(U n ) is nonsingular, 2 and positive definite 3 , [36][37][38]. This means that B(U n ) is nonsingular and negative definite, i.e., all the eigenvalues are negative, and then λ(I − t x 2 B(U n )) > 1. This confirms the statement. Although this analysis is only performed on first order SI scheme, numerical experiments in Sect. 3 demonstrate its validity for higher-order SI-IMEX-RK schemes as well.

High Order Semi-implicit Temporal Discretization
We introduce our proposed semi-implicit strategy for the time discretization in the framework of high order IMEX-RK methods. We emphasize that the following strategy avoid of solving nonlinear equations required by a fully implicit scheme and avoid stringent time stepping constraint required by an explicit one. [36,37]. 2 A diagonally dominant symmetric matrix A with nonnegative diagonal entries is positive semidefinite. Moreover, if A is irreducible, and |a ii | > n j=1, j =i |a i j | for at least one i, then A is nonsingular, [37]. 3 A n × n matrix A is positive define iff it is a nonsingular positive semidefinite matrix, [37].
In the following section, we will generalize the basic idea presented in the introduction to Eq. (12). This will be achieved through the use of the approach outlined in [23] which is somehow inspired by partitioned Runge-Kutta methods, [39].
We start to consider the general class of autonomous problems of the form where F : R n × R n → R n is a sufficiently regular mapping. We observe that we can rewrite system (16) as a partitioned one with initial conditions u * (t) = u(t) = u 0 . In such case the solution of system (17) satisfies for any t ≥ t 0 and is also a solution of Eq. (16). System (17) is a particular case of partitioned system [39], with an additional computational cost since we double the number of variables. Let us remark that the duplication of the unknowns in (17) does not take place if appropriate choices of the IMEX-RK scheme are considered, [23]. Now we are ready to introduce the semi-implicit strategy for solving autonomous problems of the form (17). This strategy involves treating the first variable u * explicitly and the second variable u implicitly. From now on we shall adopt IMEX R-K schemes and the coefficients of the method are usually represented in a double Butcher tableau as For the implicit part, diagonally implicit R-K (DIRK) schemes are often employed. This approach is simpler to implement and guarantees that the terms involving the first argument of F are explicitly computed.
Then the SI-IMEX-RK method applied to (17) is implemented as follows. First we set u * n = u n and compute the stage values for i = 1, . . . , s, with the numerical solution where From now on we shall adopt IMEX R-K schemes withb i = b i , for i = 1, . . . , s. We observe that becauseb i = b i , for i = 1, . . . , s then the numerical solutions are the same, i.e. if u * 0 = u 0 then u * n = u n , for all n ≥ 0, therefore the duplication of the system is only apparent and not necessary if we adopt the RK fluxes K i = F(U * i , U i ) as basic unknowns, so that one can rewrite the scheme in the following form for i = 1, . . . , s, with the numerical solution In light of the previous discussion, the main advantage of using the SI approach to compute numerical solutions is to solve a linear system. In the case of a high-order SI-IMEX-RK scheme, the system (21) is linear in K i . For instance, the system (12) can be expressed as linear equations in terms of K i : are computed explicitly. Now, we provide the expression for the function F(u * , u) that applies to Eqs. (1), (2), (6), (7) and (8). Then, we have for the diffusion equations, whereas for the dispersion equation (6), rewriting the nonlinear Finally, for the nonlinear term in the fourth order diffusion equations (7) and (8) we have So in practice, after discretizing the spatial operators of Eqs. (23), (24), and (25) on a chosen space grid with they are converted into a semi-discrete form. Specifically, for Eqs. (1)- (7), we have Then the resulting SI-IMEX-RK scheme can be applied. We conclude this section by noting that it will be useful from now on to characterize different IMEX schemes according to the structure of the DIRK method. Following [40] we say an IMEX-RK method is of type I if the matrix A ∈ R s×s is invertible, and we say an IMEX-RK method is of type II [14] if the matrix A can be written as with a = (a 21 , . . . , a s1 ) T ∈ R (s−1) and the submatrixÂ ∈ R (s−1)×(s−1) is invertible, or equivalently a ii = 0, i = 2, . . . , s and the DIRK method is reducible to a method using s − 1 stages. In the special case a = 0, b 1 = 0 the scheme is said to be of type ARS [13]. Finally, we list below the third order IMEX-RK schemes that we shall use in the paper. The triplet (s, σ, p) characterizes the number of stages of the implicit scheme s, the number of stages of the explicit scheme σ and the order of the scheme p.
where α = 0.24169426078821, β = α/4 and η = 0.12915286960590. We call this scheme SSP-DIRK3(4,3,3). This scheme is of type I. Some coefficients of the explicit part are computed by assuming that the corresponding stability region is the same of a classical fourth stage, fourth-order explicit RK schemes, i.e., The stability functions (31) guarantees to have a larger stability region for the explicit part of the scheme than the one of SSP-DIRK3(4,3,3) scheme (29).
with initial conditions u * (t) = u(t) = u 0 , and applying a SI-IMEX-RK scheme we have with the numerical solution Generally, K i and L i given by Eq. (34) for 1 ≤ i ≤ s are different. However, there are two cases where K i = L i for i = 1, . . . , s. The first case is when the system is autonomous, which means that F does not explicitly depend on time, i.e., Eq. (17). The second case is whenc i = c i for i = 1, . . . , s. In these two cases, only one evaluation of F is needed in Eq. (34) and only one set of stage fluxes is computed, i.e., (21).
Note that this restriction can be removed even in the general case of a non-autonomous system, andc i = c i , for i = 1, . . . , s, it is still possible to derive a scheme that does not require two sets of stage fluxes (34). The details are reported in the "Appendix" of the paper [23]. Now we observe that selecting a different vector of weights for the U * , sayb i = b i , for i = 1, . . . , s, such a vector will provide a lower/higher order approximation of the solution for U * . This would allow for automatic time step control. This procedure is commonly used in numerical methods for ODEs [39]. In practice, we use vectorb for evaluating the U * variable, and by setting u * n = u n at the beginning of each time step, one advances the numerical solution with the more accurate one U , and uses the other variable U * to estimate the error. However, in this paper, we do not implement any time step control.
Finally looking at the list of third-order IMEX-RK schemes presented previously, we observe that schemes ARS (3,4,3) and ARS(4,4,3) havec i = c i , and therefore require only one evaluation of F, i.e., (21). On the other hand, schemes SSP-DIRK3(4,3,3) and I-IMEX (3,4,3) have c 1 = 0, andc i = c i , for all i = 2, . . . , s. Consequently, K 1 = L 1 for the first evaluation of F whereas K i = L i for i = 2, . . . , s. Nonetheless, these two schemes haveb i = b i for all i = 2, . . . s andb 1 = b 1 = 0, which implies that the numerical solutions for both variables U * and U are identical because we do not need to use the stage fluxes K 1 and L 1 in the evaluation of the numerical solution. Additionally, scheme ARS(4,4,3) has b i =b i for all i. To handle this case, we use the variable U to compute the numerical solution. Specifically, at the beginning of each time step n, we set u * n = u n , and then solve using (21), since c i =c i for all i, and then (22) for the numerical solution.

The Spatial Discretization
Some preliminary notations are given. We consider one dimensional problems and assume that the computational domain [a, b] is uniformly partitioned into M cells, with spatial mesh size x = (b − a)/M. We use the lowercase letter u i to denote the numerical solution to the equation at the the spatial position x i for i = 1, . . . , M.

The Spatial Discretization for the Diffusion Equation
For the spatial discretization of the nonlinear diffusion equation we take the following formula that provides a fourth order approximation to (a(u)u x ) x , with a five-point stencil [41,42] Note that in the case a(u) = 1 the formula (35) becomes the classical fourth order central finite difference scheme, i.e., D 2

The Spatial Discretization for the Dispersive Equation
To discuss about high order space discretization for Eq. (6), and in particular, for the nonlinear term (u(u n x x )) x , we consider again the equivalent nonlinear term: with a(u) = nu n−1 . Then, to approximate the term (37) we use the following fourth order finite difference spacial approximation, where is a fourth order central approximation for u x (x j ).

The Spatial Discretization for the Biharmonic-type Equation
The fourth order finite difference scheme for the fourth order spatial derivative in the biharmonic-type equation (7) can be written as where D 2 x u j is a fourth order centered difference approximation to u x x (x j ) defined as (36). Note that in (7) when a(u x ) = 1 we get from (39), D 2 x D 2 x u j , i.e., the central difference approximation for the fourth order derivative u x x x x :

Numerical Results
In this section, we will numerically verify the orders of accuracy and performance of SI-IMEX-RK schemes for the solution of time dependent PDEs with high order derivatives, as Eqs. (1), (2), (6), (7) and (8). For simplicity we used in all of our examples periodic boundary conditions, although most of our discussions can be adapted for other types of boundary conditions. In each of the numerical tests presented below, we consider the third order IMEX-RK schemes introduced in Sect. 2.2 for the time discretization and we provide the time step t that will produce stable solutions across the entire solution domain. For the space discretization we use high-order finite difference schemes introduced in Sect. 2.3, and we discretize the convection term by means of standard third order finite difference WENO schemes with local Lax-Friedrichs flux, [43][44][45]. All the tests presented in this paper are performed in one-dimension since the main issue of this paper is to design an appropriate high order time discretization for time-dependent PDEs. Finite difference schemes have been extensively applied for solving higher dimensional equations because of its simplicity in design and coding [46][47][48], then it is straightforward to extend to higher-dimensional equations this type of space discretization.
It is important to note here that, as mentioned in Remark 2.1, although some equations in the test examples below are presented in non-autonomous form due to the time-dependence of the source term, the errors tables reported below are computed solely based on the numerical solution obtained through (22) using the SI framework (21)- (22), where the coefficients are given by the IMEX-RK tables in Sect. 2.2.

The Second Order Diffusion Equations
Test 1. First we consider the diffusion equation Here we consider the case a(u) = u 2 + 1 with initial condition u(x, 0) = sin(x) and the source term f (x, t) is chosen such that the exact solution is u(x, t) = sin(x − t).
For this example, we take: the time step t = C x, with C = 1 and x = 2π/N , periodic boundary conditions and final time T = 10. In Table 1 we display the numerical results of the third order IMEX-RK schemes given in Sect. 2.2 coupled with the finite difference type spatial discretization (35). In the table we list the L 2 , L 1 and L ∞ errors and the orders of accuracy of each scheme and we can see that the schemes are stable and achieve the third order of accuracy. Now we use a larger value of C to show that the schemes remain stable even for values of C > 1. In the Table 2 special attention has been paid to large C = 10, where the numerical results can be compared to those of C = 1. As example, we report the results for the schemes: I-IMEX(3,4,3) (30) and ARS(3,4,3) (27), the other schemes produce similar results. We see that the schemes are stable in both cases and achieve optimal orders of accuracy. Note that larger C, as in the case C = 10, causes larger errors. Test 2. Next we consider the convection-diffusion equation where the diffusion coefficient is a(u) = u 2 + 2, initial conditon u(x, 0) = sin(x) and the source term f (x, t) = 1 4 (4 cos(x + t) + 9 sin(x, t) + 2 sin(2(x + t)) − 3 sin(3(x + t))) .
The problem has the exact solution: u(x, t) = sin(x + t).
The standard third order WENO(3,2), [45] scheme with linear weights is used for the discretization of the convection term (u 2 /2) x . We compute to T = 1 with the time step t = x. The numerical results with different IMEX-RK schemes are listed in Table 3. From the experiment we can see that SI strategy based on IMEX-RK schemes coupled with finite difference schemes are stable and achieve optimal orders of accuracy. Test 3. Now we consider the porous media equation (PME) in which m is a constant greater than one. This equation describes a lot of diffusion processes, such as in nonlinear problems of heat and mass transfer, combustion theory, and flow in porous media. where u is either a concentration or a temperature required to be non negative. The initial condition u 0 (x) is assumed to be bounded non negative continuous function. Equation (42) can be written as with a(u) = mu m−1 , [26]. It's a degenerate parabolic equation and it degenerates at points where u = 0, resulting in the phenomenon of finite speed of propagation and sharp fronts.
Here, we consider PME (43) with the Barenblatt solution (44). As initial condition we take the (44) at time t = 1 and we begin the computation from t = 1 in order to avoid the singularity near t = 0. The boundary condition is u(±6, t) = 0 for t > 1 and we use uniform mesh with N = 160 points and time step t = x. We adopt for the spatial discretization the fourth order approximation (35), and for the time discretization the third-order I-IMEX(3,4,3) scheme (30). In Fig. 1we plot the numerical results for m = 2, 3, 5 at t = 2 and we see that the numerical solution simulates the Barenblatt solution (44) accurately and sharply. However, there is a slight undershoot for the numerical solution for m = 3, 5 near the discontinuity. We verify that this undershoot is reduced when we take more mesh points, for example N = 320 (see Fig. 1d). We point out that in this simulation we adopted a not suitable space discretization for the PME equation (42), in the literature several optimal discretization are introduced for PME equation, for instance, in [58] an adaptation of the WENO technique has been proposed in oder to maintain conservation, accuracy and non-oscillatory performance. Similarly in [26] the authors used local discontinuous Galerkin methods coupled with two specific explicit-implicit-null time discretization for solving onedimensional nonlinear diffusion problems (43). However, the I-IMEX(3,4,3) scheme (30) produces numerical solutions accurately without noticeable oscillations for the PME equation (42). From the above test, we can conclude that the SI-IMEX-RK approach essentially: (1) removes the strict parabolic time step restriction, i.e. t ≈ x 2 , that usually one has using explicit time discretization methods for nonlinear diffusion problems (43) and (2) avoids nonlinear solvers required using implicit time discretization methods. Test 4. We conclude with an example of a strongly degenerate parabolic convectiondiffusion equation presented in [59]: We take = 0.1, f (u) = u 2 and The equation is therefore of hyperbolic nature when u ∈ [−0.25, 0.25] and parabolic elsewhere. We solve this problem considering the following initial data Here we used the I-IMEX(3,4,3) scheme (30) for the integration in time and we consider the classical hyperbolic CFL condition to set t. In Fig.2 , the numerical simulations for different numbers of grid points N =
We compute to T = π, with λ = 0.1. We choose t = x for the I-IMEX(3,4,3) (30) and SSP-DIRK3(4,3,3) (29) whereas, as mentioned in the introduction about the stability, for the scheme ARS(3,4,3) (27) in order to avoid that the numerical solution blows up we reduced the time step to t = 0.5 x. In Table 4 we show the convergence rate for these schemes. We omit the convergence results of the ARS(4,4,3) scheme (28) in this table because they are comparable to those of the ARS(3,4,3) scheme. As expected, the SI-IMEX-RK finite difference schemes are stable and achieve the correct order of accuracy. Note that, since λ represents the velocity of the traveling wave, selecting a larger value of λ requires a careful choice of the time step to ensure the stability of the method. In such cases, we set the time step as with a fixed CFL. As an example, we provide Table 5 below that shows the convergence rates of the schemes I-IMEX (3,4,3) and ARS (3,4,3) with λ = 10 and C F L = 0.5. We see that the schemes are stable and achieve the expected order of accuracy. Note that larger λ causes larger errors. We do not include the convergence results of the other schemes Table 4 The   Table 5 The  since they are comparable to those presented in Table 5. In all of the following tests: 6., 7., and 8., the boundary conditions are taken to be periodic in an interval much larger than the compact support of the initial conditions. Test 6. Now, we consider the nonlinear dispersive equation (6) with m = 2, and n = 1. In this example we consider the following compacton initial data [10] u(x, 0) = 2 cos 2 (x/2) |x| ≤ π, 0, |x| > π .
In this case the numerical solution is a traveling wave given by with λ = 1. Here, in order to ensure the stability of the schemes, we consider (48) with to set t. In Table 6 we show the convergence rate of the schemes: ARS (3,4,3), SSP-DIRK(4,3,3) and I-IMEX (3,4,3) and is approximately three for L 1 , L 2 , and L ∞ -norm. Again, we omit the convergence results of the ARS(4,4,3) scheme in this table because they are comparable to those of the ARS(3,4,3) scheme. Test 7. In this example we approximate solutions of Eq. (49) by showing a collision for three compactons, [6]. The initial data are taken as The computational domain is [−4π, 22π] and N = 300 mesh points. Here we choose the third order I-IMEX(3,4,3) scheme (30) in time.
The results are shown in Fig. 3 . As one can see rapid oscillations develop behind the moving compactons and we would like to get a scheme that is at least stable with respect to small perturbations of the solution. The SI approach coupled with WENO space discretization does not solve correctly the oscillations in the tails but it is able with a large time step t (48) to control the oscillations and ensure that the solution does not blow up for a long time. Note that suitable space discretization [6,10] are capable of capturing simultaneously the oscillations and the non-smooth fronts. For example, in [6] the authors presented an approach for the construction of the numerical fluxes in which the scheme is stable with respect to small perturbations of the solution and then to avoid oscillations in the tails. This latter case is beyond the goal of this paper. Test 8. We conclude this section with the following test, [10]. We consider Eq. (6) with m = 2 and n = 1. We use arbitrary initial data In this case we take N = 200 and t given by (48). In time, we use the third order I-IMEX (3,4,3) scheme (30). In Fig. 4we plot the solution at time T = 0, 1, 2, 4, 6, 8.  As we can see in the figures the behaviour is similar to the results reported in [10]. Clearly, the results of the SI approach coupled with WENO discretization suffer from spurious oscillations in the tail for T = 4, 6, 8. Even in this case, our approach does not solve correctly the oscillations in the tail but it is able with a large time step t (48) to control the oscillations and ensure that the solution does not blow up for a long time. In order to avoid the numerical oscillations in the tail, suitable space discretization can be used as in [10] where a particle method has been developed.

Fourth Order Diffusion Equation
Test 9. First we consider the following nonlinear biharmonic-type equation [20] u t + (( with initial condition u(x, 0) = sin(x), source term and exact solution u(x, t) = e −t sin(x). We compute the solution to T = 1 with the time step t = C x setting C = 1 and C = 6.5. The numerical errors and orders of accuracy are listed in Tables 7and 8. In each table, we display the numerical results for the schemes SSP-DIRK(4,3,3) (29) and I-IMEX(3,4,3) (30). From the experiment we can see that the schemes are stable and achieve optimal orders of accuracy. However, in the two tables the I-IMEX(3,4,3) scheme gives a better orders of accuracy than SSP-DIRK(4,3,3) scheme. This verifies that the better condition on the stability (see Sect. 2.2) for the scheme I-IMEX(3,4,3) guarantees a good behaviour with mesh refinements.
We also note that in the Table 8 even if larger C causes larger error than C = 1, the order of accuracy of the schemes achieves the correct order of convergence. Test 10. Next, we have the biharmonic-type equation with a flux term: We take initial condition u(x, 0) = sin(x) and the source term f (x, t) is chosen such that the exact solution is We compute the solution to T = 1 with the time step t = x. The numerical errors and orders of accuracy are listed in  capture correctly the order of accuracy. We omitted these results here to save space. We remember that generally for the 4-th order PDEs, an explicit RK method requires severe time step restriction, i.e., t = O( x 4 ) for stability. In the future, we aim to conduct a more detailed stability analysis of these ARS schemes, which would require additional details. It is worth noting that an interesting aspect of this analysis could be the structure of the matrix in the the implicit part of the ARS scheme, i.e., a 11 = 0 in (26), which could potentially result in instability in the subsequent internal stages of the explicit part after the first stage. However, this is currently only a tentative conjecture, and we will postpone its investigation to a future work.

Concluding Remarks
In this paper, following the idea of the work [23] we have developed a SI strategy based on IMEX-RK methods coupled with high order finite difference schemes for solving high order dissipative, dispersive and special biharmonic-type equations in one dimension. The SI IMEX-RK schemes so designed for high order time-dependent PDEs does not need any nonlinear iterative solver, and not require any time step restriction that usually one has using explicit methods. Numerical experiments show that the schemes are stable and achieve the aspected orders of accuracy for large time step.
Since the main issue of this paper is to design an appropriate high-order time discretization for time-dependent PDEs, the numerical experiments are only performed for one-dimensional equations. For the space discretization, we considered only classical finite difference spatial discretization because of its simplicity in design and coding and it is straightforward to extend to higher-dimensional equations. However, other types of space discretization can be used.

Author Contributions Formal analysis, Methodology, Software
Funding Open access funding provided by Università degli Studi di Catania within the CRUI-CARE Agreement. Sebastiano Boscarino is partially supported by the Italian Ministry of Instruction, University and Research (MIUR), with funds coming from PRIN Project 2017 (2017KKJP4X, entitled "Innovative numerical methods for evolutionary partial differential equations and applications"). The authors is a member of the INdAM Research group GNCS. Furthermore, Sebastiano Boscarino is supported for this work by the Spoke 1 "FutureHPC & BigData" of the Italian Research Center on High-Performance Computing, Big Data and Quantum Computing (ICSC) funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di "campioni nazionali di R&S (M4C2-19 )".

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The author declares that he has no known competing financial interests that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.