Multiscale RBF-based central high resolution schemes for simulation of generalized thermoelasticity problems

Abstract

In this study, average-interpolating radial basis functions (RBFs) are successfully integrated with central high-resolution schemes to achieve a higher-order central method. This proposed method is used for simulation of generalized coupled thermoelasticity problems including shock (singular) waves in their solutions. The thermoelasticity problems include the LS (systems with one relaxation parameter) and GN (systems without energy dissipation) theories with constant and variable coefficients. In the central high resolution formulation, RBFs lead to a reconstruction with the optimum recovery with minimized roughness on each cell: this is essential for oscillation-free reconstructions. To guarantee monotonic reconstructions at cell-edges, the nonlinear scaling limiters are used. Such reconstructions, finally, lead to the total variation bounded (TVB) feature. As RBFs work satisfactory on non-uniform cells/grids, the proposed central scheme can handle adapted cells/grids. To have cost effective and accurate simulations, the multiresolution–based grid adaptation approach is then integrated with the RBF-based central scheme. Effects of condition numbers of RBFs, computational complexity and cost of the proposed scheme are studied. Finally, different 1-D coupled thermoelasticity benchmarks are presented. There, performance of the adaptive RBF-based formulation is compared with that of the adaptive Kurganov-Tadmor (KT) second-order central high-resolution scheme with the total variation diminishing (TVD) property.

Graphical abstract

Keywords

central high resolution schemes / RBFs / higher order accuracy / generalized thermoelasticity / multiresolution-based adaptation

1 Introduction

A wide range of applications of thermoelastic behaviors in applied engineering (such as pressure vessels, pipes and aerospace), engineering science (e.g., pulse laser irradiations and rapid solidifications) and natural phenomena (for example seismology) lead to extensive attention on this research area in recent years [1,2]. Regarding the thermoelastic phenomenon, a good survey of some applications and relating numerical methods can be found in [3]. Interaction of thermo-mechanical responses can considerably affect structural responses and even may lead to the failure. These behaviors are important especially for heat shocks due to sudden and considerable arising of stresses in solids. Some experiments reveal that heat shocks propagate with finite speeds in an elastic medium when the heat conduction takes place at a short time interval [4]. This is in contrast with classical thermoelastic theories predicting an infinite propagating speed. To explain this phenomenon, different theories have been developed. Lord and Shulman (LS) [5] introduced the theory of generalized thermoelasticity with one relaxation time for isotropic bodies, based on the Cattaneo heat conduction law. Green and Lindsay [6] developed the generalized thermoelasticity with two relaxation times. Green and Naghdi (GN) [7] proposed the theory of thermoelasticity without energy dissipation. They also developed two other theories [8,9]; for a general review also see [10]. These coupled formulations were also extended to consider fractional equations of the generalized thermoelasticity theories [11,12]. For practical problems, later, generalized thermoelasticity theories with variable material properties have been investigated [13–26].

Aforementioned generalized thermoelastic theories lead to hyperbolic systems where shock waves (singular propagating fronts) develop in their solutions. Regarding mechanical waves, discontinuous (or singular) solutions can be developed due to: 1) material nonlinearities; 2) abrupt changing in material properties, such as: existence of a narrow fluid-filled crack [27], random fluid-filled cavities [28], stratified media [27] or composite materials [22,24,25,29‒31]; 3) loading properties; e.g., those used for simulation of seismic sources [32] or stress shocks deduced by abrupt thermal loadings; 4) equation nonlinearity due to coupling with thermal or magneto-thermal effects [10]. In most of numerical methods, existence of such discontinuous fronts can lead to numerical (artificial) dissipation-dispersion and finally numerical instabilities. High-resolution schemes have successfully been developed for numerical simulations of hyperbolic problems. Godunov [33] showed that linear methods cannot provide non-oscillatory solutions higher than one. In this regard, different nonlinear methods have been developed; some of the most famous approaches are: 1) essentially non-oscillatory (ENO) [34], and weighted-ENO (WENO) schemes [35], 2) wave-propagation algorithms [36]; 3) central schemes [37,38] (based on the MUSCL methods using the concept of the flux/slope limiter [39]); and 4) central-WENO schemes [40–42]. Based on total variation (TV) of these schemes, they are mostly total variation diminished (TVD) or total variation bounded (TVB). These criteria are to guarantee numerical stability; they prevent or control developing of spurious oscillations in numerical solutions. TVD-based schemes have two limitations: 1) they mostly have second order accuracy on smooth solutions and first order accuracy around discontinuities, 2) designing of higher-order TVD schemes are restrictive. To cure these shortcomings, TVB-based methods have been developed. Therefore, this criterion is used in this study to develop a stable adaptive higher-order scheme.

High resolution schemes have been used for numerical simulations of generalized thermoelastic equations. Levy et al. used a TVD solver for such implementation in porous media and compared numerical and experimental results [43,44]. Another TVD-based algorithm was used for simulation of elastic-plastic solids with the thermoelastic behavior [45]. The wave-propagation algorithm was also utilized for simulation of thermoelastic waves [46,47]. It should be mentioned that above-mentioned studies were performed on uniform cells/grids.

In this study, the central high resolution schemes are used, as they offer the following benefits: 1) identifying of propagation directions is not necessary. This is a complex procedure, especially for 2-D and 3-D problems; 2) they are Riemann-free solvers; 3) both staggered and non-staggered solvers were developed; 4) both fully-discrete and semi-discrete forms were provided; 5) both TVD and TVB based formulations were developed; 6) second and higher order formulations were developed. In summary, the central high resolution schemes can act as a black-box for numerical simulations of nonlinear first-order hyperbolic equations.

Most of the high resolution schemes use polynomials for reconstruction of state variables over cells. Other candidates with desired features have also been proposed. The family of radial basis functions (RBFs) is one of them; they have successfully been used for several problems with advection-diffusion properties, e.g., Refs. [48‒52]. Their optimum recovery feature (from given cell averages) was proved in [53], where the concept of the optimum recovery is defined in [54,55]. RBF-based interpolations also lead to reconstructions with minimized roughness [56]. These features are important for stable and oscillation-free solvers and make RBFs as a perfect candidate for the reconstructions. RBFs have successfully been used in the ENO [50,53,57] and WENO [50,52,58,59] high resolution schemes of higher order accuracy.

In this regard, in this study, RBFs are integrated with the central high resolution schemes to provide stable higher order solvers. Average-interpolating RBFs are used for reconstruction of state variables (over each cell). To guarantee the monotonicity at cell edges, nonlinear scaling limiters are utilized. Using of these limiters along with smoothness feature of RBFs lead to TVB central high resolution solvers. As RBFs can effectively handle non-uniform cells/grids, corresponding central high resolution schemes can properly work on adaptive cells/grids. And this is an important capability, as high-order high resolution schemes are expensive. On the other hand, by the grid adaptations, condition numbers in RBF formulations can increase considerably around finer grids; this leads to growing of errors and even numerical instabilities. Therefore, there is a trade-off between accuracy of adapted solutions and corresponding errors [60]. Shape parameters of RBFs control performance of these functions. So at first a general survey is presented for different selection methods of shape parameters and then a simple approach is followed for such choosing.

Here, cells/grids are adapted by the multiresolution analysis (MRA). Coefficients of wavelet transforms have two properties [61]:

1) They concentrate automatically around high gradient zones and discontinuities; this is a perfect criterion for capturing of such zones,

2) There is a one-to-one correspondence between grid points and coefficients of the transform (with the pyramid algorithm).

These features make the MRA as a perfect tool for grid adaptations. For the wavelet transform, the Dubuc-Deslauriers (D-D) interpolating wavelets [61,62] are used which lead to a simple and fast transform algorithm. All calculations are directly performed in the physical domain. So, the transform coefficients have physical meanings.

Finally, it should be mentioned that multiresolution-based adaptation procedures have successfully been integrated by upwind-based high resolution schemes for simulation of hyperbolic partial differential equations (PDEs) [62–69].

This paper is composed of eight parts, Section 2 presents the unified theory of the generalized thermoelasticity with constant and variable coefficients; there, corresponding hyperbolic systems with their Jacobian matrices and eigenvalues are derived. Section 3 devotes to the main concept of MRA, MRA-based grid adaptation and MRA-based modification of adapted grids. Section 4 explains average interpolating RBFs and corresponding reconstructions. In this section, the concept of the nonlinear scaling limiters is used to monotonize RBF-based reconstructions at cell-edges to have monotonic piecewise reconstructions. Afterwards, choosing a good value for shape parameters of RBFs is discussed. Section 5 presents relationships between polynomial- and RBF-based reconstructions; based on this, accuracy order of RBF-based formulations is evaluated. Then, the semi-discrete form of such central high resolution schemes is provided. At the end of this section, computational complexity is studied and then computational cost of the proposed scheme is compared with two other third order central high resolution schemes. Section 6 describes a numerical algorithm for implementation of the RBF-based central scheme on MRA-based adapted grids. Some benchmarks are simulated in Section 7. The conclusion is presented at the end of the paper.

2 Generalized thermoelasticity and corresponding hyperbolic systems

In this section, at first, the unified-form of generalized thermoelasticity is presented. These equations are then rewritten for two 1-D cases: thermoelastic problems with constant and thermal-dependent coefficients. For each one of these cases, corresponding first-order hyperbolic system, Jacobian matrix and eigenvalues are presented. For each formulation, the maximum value of eigenvalues (regardless of propagation directions) is needed for numerical discretization with high-order central high resolution schemes.

2.1 Uniform formulation of generalized thermoelasticity

Regarding the GN and LS generalized thermoelasticity formulations for isotropic homogeneous materials, a unified representation can be expressed. For this purpose, conservation/balance laws of momentum and energy along with the constitutive equation of materials and the generalized heat conduction equitation should be presented. These relationships can be expressed as [24]:

1) The conservation of momentum:

(1){{\displaystyle \rho u}}_{i,tt}={{\displaystyle \rho f}}_i+{{\displaystyle \sigma }}_{ij,j},\hspace{0.17em}\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}i,j=\{\mathrm{1},\mathrm{2},\mathrm{3}\},

2) The conservation of energy:

(2){{\displaystyle q}}_{i,i}=\rho r-\rho {{\displaystyle c}}_p{{\displaystyle \theta }}_{,t}-{{\displaystyle T}}_{\mathrm{0}}\beta {{\displaystyle \gamma }}_{kk,t,}\hspace{0.17em}\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}i,k=\{\mathrm{1},\mathrm{2},\mathrm{3}\},

3) The constitutive equation:

(3){{\displaystyle \sigma }}_{ij}=\lambda {{\displaystyle \gamma }}_{kk}{{\displaystyle \delta }}_{ij}+\mathrm{2}{{\displaystyle \mu \gamma }}_{ij}-\beta \theta {{\displaystyle \delta }}_{ij},\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}i,j,k=\{\mathrm{1},\mathrm{2},\mathrm{3}\},

with the linear strain-displacement relationship {{\displaystyle \gamma }}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({{\displaystyle u}}_{i,j}+{{\displaystyle u}}_{j,i}),

4) The generalized heat conduction equation:

(4){{\displaystyle \eta }}_{\mathrm{1}}{{\displaystyle q}}_i+{{\displaystyle \tau }}_{\mathrm{0}}{{\displaystyle q}}_{i,t}+{{\displaystyle \eta }}_{\mathrm{2}}{{\displaystyle q}}_{i,t}=-{{\displaystyle \eta }}_{\mathrm{1}}{{\displaystyle k\theta }}_{,i}-{{\displaystyle \eta }}_{\mathrm{2}}{k}^{*}{{\displaystyle \theta }}_{,i},\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}i=\{\mathrm{1},\mathrm{2},\mathrm{3}\},

where:{{\displaystyle u}}_i, {{\displaystyle f}}_i, and {{\displaystyle q}}_i denote respectively components of displacements, body force, and heat flux vectors; {{\displaystyle \sigma }}_{ij} and {{\displaystyle \gamma }}_{ij} are elements of the stress and strain tensors; {{\displaystyle Z}}_{,i}:=\partial Z/{{\displaystyle \partial x}}_i and {{\displaystyle Z}}_{,t}:=\partial Z/\partial t; \theta =T-{{\displaystyle T}}_{\mathrm{0}} shows temperature increment, in which T and {{\displaystyle T}}_{\mathrm{0}} are the absolute and reference temperatures, respectively; \rho is the mass density; {{\displaystyle c}}_p is the specific heat at a constant strain; r is the internal heat source; \beta =(\mathrm{3}\lambda +\mathrm{2}\mu ){{\displaystyle \alpha }}_T is the thermal-mechanical coupling coefficient wherein {{\displaystyle \alpha }}_T is the coefficient of linear thermal expansion and \lambda & \mu are the Lame’s constants; k denotes the thermal conductivity; k^{*} is a new material constant introduced in the GN formulation; and {{\displaystyle \eta }}_{\mathrm{1}}, {{\displaystyle \eta }}_{\mathrm{2}} and {{\displaystyle \tau }}_{\mathrm{0}} are theory-dependent constants. For different theories, these parameters are defined as:

1) The LS model: {{\displaystyle \eta }}_{\mathrm{1}}=\mathrm{1}, {{\displaystyle \eta }}_{\mathrm{2}}=\mathrm{0} and {{\displaystyle \tau }}_{\mathrm{0}}\ne \mathrm{0},

2) The GN model:{{\displaystyle \eta }}_{\mathrm{1}}=\mathrm{0}, {{\displaystyle \eta }}_{\mathrm{2}}=\mathrm{1} and {{\displaystyle \tau }}_{\mathrm{0}}=\mathrm{0},

3) The classical coupled model:{{\displaystyle \eta }}_{\mathrm{1}}=\mathrm{1}, {{\displaystyle \eta }}_{\mathrm{2}}=\mathrm{0}, and {{\displaystyle \tau }}_{\mathrm{0}}=\mathrm{0}.

Introducing the velocity v={{\displaystyle u}}_{,t} as an independent variable, Eq. (1) can be rewritten as the velocity-stress equations:

(5){{\displaystyle \rho v}}_{i,t}={{\displaystyle \rho f}}_i+{{\displaystyle \sigma }}_{ij,j},

where this system is the first-order form of Eq. (1).

One dimensional thermoelasticity with constant coefficients

In the following, it is assumed that there is no variations in u, v, \sigma, {{\displaystyle \theta }}_{,} and q in the y and z directions: an 1-D problem, where i=j=\mathrm{1}. With this assumption, formulations of corresponding system of 1-D hyperbolic equations are provided. For this, at first, let us calculate {{\displaystyle \theta }}_{,t}:={{\displaystyle \theta }}_{\mathrm{1},t} from Eq. (2) and also take the first derivative of Eq. (3) respect to the time; in this new equation substitute {{\displaystyle \theta }}_{,t} (obtained from Eq. (2)) to attain:

(6){{\displaystyle \sigma }}_{,t}=(\lambda +\mathrm{2}\mu )(\mathrm{1}+v){{\displaystyle v}}_{,x}+{\displaystyle \frac{\beta }{{{\displaystyle \rho c}}_p}}q,x-{\displaystyle \frac{\beta }{{{\displaystyle c}}_p}},

where v:={{\displaystyle v}}_{\mathrm{1}}, \sigma :={{\displaystyle \sigma }}_{\mathrm{1}}, q:={{\displaystyle q}}_{\mathrm{1}} and v:={\displaystyle \frac{{{\displaystyle T}}_{\mathrm{0}}{\beta }^{\mathrm{2}}}{{{\displaystyle \rho c}}_p[\lambda +\mathrm{2}\mu ]}} denotes the thermal coupling constant.

Eqs. (5), (6), (2), and (4) can be rewritten as:

(7){{\displaystyle \mathit{U}}}_{,t}+{{\displaystyle \mathit{F}}}_{,X}=\mathit{S},

with the semi-linear form:

(8){{\displaystyle \mathit{U}}}_{,t}+\mathit{A}.{{\displaystyle \mathit{U}}}_{,x}=\mathit{S},

where: U, F and S are vectors of state variables, fluxes and sources, respectively; A:={\displaystyle \frac{\partial A}{\partial U}} (with components {{\displaystyle A}}_{ij}:={\displaystyle \frac{{{\displaystyle \partial A}}_i}{{{\displaystyle \partial U}}_j}}) denotes the Jacobian matrix. For the state vector (U) and the flux vector (F)

\mathit{U}={\{v,\sigma ,q,\theta \}}^T,

(9)\mathit{F}={{\displaystyle \left\{-{\displaystyle \frac{\mathrm{1}}{\rho }}\sigma ,-{{\displaystyle v}}_e^{\mathrm{2}}\rho [\mathrm{1}+v]v-{\displaystyle \frac{\beta }{{{\displaystyle \rho c}}_p}}q,\overline{a}\theta ,{\displaystyle \frac{v\rho {{\displaystyle v}}_e^{\mathrm{2}}}{\beta }}v+{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \rho c}}_p}}q\right\}}}^T,

the Jacobian matrix becomes:

(10)\mathit{A}=\left(\begin{array}{llll}\mathrm{0}& -{\displaystyle \frac{\mathrm{1}}{\rho }}& \mathrm{0}& \mathrm{0}\\ -{{\displaystyle v}}_e^{\mathrm{2}}\rho [\mathrm{1}+v]& \mathrm{0}& -{\displaystyle \frac{\beta }{{{\displaystyle \rho c}}_p}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \overline{\alpha }\\ {\displaystyle \frac{{{\displaystyle v\rho v}}_e^{\mathrm{2}}}{}}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{{{\displaystyle \rho c}}_p}}& \mathrm{0}\end{array}\right),

and the source vector is:

(11)S={{\displaystyle \left\{f,-{\displaystyle \frac{\beta }{{{\displaystyle c}}_p}}r,-{\displaystyle \frac{{{\displaystyle \eta }}_{\mathrm{1}}}{{{\displaystyle \eta }}_{\mathrm{2}}+{{\displaystyle \tau }}_{\mathrm{0}}}}q,{\displaystyle \frac{\mathrm{1}}{{{\displaystyle c}}_p}}r\right\}}}^T,

where {{\displaystyle v}}_e:=\sqrt{{\displaystyle \frac{\lambda +\mathrm{2}\mu }{\rho }}} is the velocity of uncoupled elastic wave and \bar{\alpha }:={\displaystyle \frac{{{\displaystyle \eta }}_{\mathrm{1}}k+{{\displaystyle \eta }}_{\mathrm{2}}{k}^{*}}{{{\displaystyle \eta }}_{\mathrm{2}}+{{\displaystyle \tau }}_{\mathrm{0}}}}. Eigenvalues of A are:

(12){{\displaystyle \bar{\lambda }}}_{\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}}=\pm {\displaystyle \frac{{{\displaystyle v}}_e}{\sqrt{\mathrm{2}}}}\sqrt{\left({\displaystyle \frac{{{\displaystyle \bar{\alpha }v}}_t^{\mathrm{2}}}{k^{*}}}+v+\mathrm{1}\right)\pm \sqrt{{\displaystyle \frac{{{\displaystyle \mathrm{4}v\bar{\alpha }v}}_t^{\mathrm{2}}}{k^{*}}}+{{\displaystyle \left({\displaystyle \frac{{{\displaystyle \bar{\alpha }v}}_t^{\mathrm{2}}}{k^{*}}}-v-\mathrm{1}\right)}}^{\mathrm{2}}}},

where {{\displaystyle v}}_t:=\sqrt{{\displaystyle \frac{\lambda +\mathrm{2}\mu }{\rho }}} denotes the velocity of uncoupled thermal wave.

One dimensional thermoelasticity with variable coefficients

In this subsection, it is assumed that the thermal parameters (k, k^{*} and {{\displaystyle c}}_p) depend linearly to the temperature increment \theta, as:

k:=k(\theta )={{\displaystyle k}}_{\mathrm{0}}(\mathrm{1}+{{\displaystyle \chi }}_{\mathrm{1}}\theta ),

(13)k^{*}:=k^{*}(\theta )={{\displaystyle k}}_{\mathrm{0}}^{*}(\mathrm{1}+{{\displaystyle \chi }}_{\mathrm{2}}\theta ),

{{\displaystyle c}}_p:={{\displaystyle c}}_p(\theta )={{\displaystyle c}}_{p\mathrm{0}}(\mathrm{1}+{{\displaystyle \chi }}_{\mathrm{1}}\theta ),

where: {{\displaystyle k}}_{\mathrm{0}}, {{\displaystyle k}}_{\mathrm{0}}^{*} and {{\displaystyle c}}_{p\mathrm{0}} are determined at the reference temperature {{\displaystyle T}}_{\mathrm{0}}; {{\displaystyle \chi }}_{\mathrm{1}} and {{\displaystyle \chi }}_{\mathrm{2}} denote small constants.

Substituting Eq. (13) into Eq. (2) yields:

(14){{\displaystyle q}}_{,x}+{{\displaystyle \rho c}}_{p\mathrm{0}}(\mathrm{1}+{{\displaystyle \chi }}_{\mathrm{1}}\theta ){{\displaystyle \theta }}_{,t}+{{\displaystyle T}}_{\mathrm{0}}{{\displaystyle \beta v}}_{,x}=\rho r,

or:

(15){{\displaystyle q}}_{,x}+{{\displaystyle \rho c}}_{p\mathrm{0}}{{\displaystyle \left(\theta +{\displaystyle \frac{{{\displaystyle \chi }}_{\mathrm{1}}{\theta }^{\mathrm{2}}}{\mathrm{2}}}\right)}}_{,t}+{{\displaystyle T}}_{\mathrm{0}}{{\displaystyle \beta v}}_{,x}=\rho r.

This equation can be written as:

(16){{\displaystyle \Theta }}_{,t}+{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \rho c}}_{p\mathrm{0}}}}{{\displaystyle q}}_{,x}+{\displaystyle \frac{({{\displaystyle T}}_{\mathrm{0}}\beta )}{{{\displaystyle \rho c}}_{p\mathrm{0}}}}{{\displaystyle v}}_{,x}={\displaystyle \frac{r}{{{\displaystyle c}}_{p\mathrm{0}}}},

where \Theta ：=\theta +{\displaystyle \frac{{{\displaystyle \chi }}_{\mathrm{1}}{\theta }^{\mathrm{2}}}{\mathrm{2}}} and so \theta =(\sqrt{\mathrm{1}+\mathrm{2}{{\displaystyle \chi }}_{\mathrm{1}}\Theta }-\mathrm{1})/{{\displaystyle \chi }}_{\mathrm{1}}.

Taking the first derivative of the constitutive equation with respect to time, and after some simplifications, we have:

(17){{\displaystyle \Sigma }}_{,t}-(\lambda +\mathrm{2}\mu ){{\displaystyle v}}_{,x}=\mathrm{0},

where \Sigma :=\sigma +\beta \theta.

Inserting Eq. (13) into the heat conduction equation (Eq. (4)), and after some simplifications, this equation becomes:

(18){{\displaystyle q}}_{,t}+{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \tau }}_{\mathrm{0}}+{{\displaystyle \eta }}_{\mathrm{2}}}}{{\displaystyle \left({{\displaystyle \eta }}_{\mathrm{1}}{{\displaystyle k}}_{\mathrm{0}}\Theta +{{\displaystyle \eta }}_{\mathrm{2}}{{\displaystyle k}}_{\mathrm{0}}^{*}\left[\Theta +{\displaystyle \frac{({{\displaystyle \chi }}_{\mathrm{2}}-{{\displaystyle \chi }}_{\mathrm{1}})(\mathrm{1}+{{\displaystyle \chi }}_{\mathrm{1}}\Theta -\sqrt{\mathrm{1}+\mathrm{2}{{\displaystyle \chi }}_{\mathrm{1}}\Theta })}{{{\displaystyle \chi }}_{\mathrm{1}}^{\mathrm{2}}}}\right]\right)}}_{,x}=-{\displaystyle \frac{{{\displaystyle \eta }}_{\mathrm{1}}}{{{\displaystyle \tau }}_{\mathrm{0}}+{{\displaystyle \eta }}_{\mathrm{2}}}}q.

The momentum equation, Eq. (1), can be written as: {{\displaystyle \rho v}}_{,t}-{{\displaystyle \sigma }}_{,x}-{{\displaystyle \beta \theta }}_{,x}+{{\displaystyle \beta \theta }}_{,x}=\rho f; since \Sigma :=\sigma +\beta \theta and \theta =(\sqrt{\mathrm{1}+\mathrm{2}{{\displaystyle \chi }}_{\mathrm{1}}\Theta }-\mathrm{1})/{{\displaystyle \chi }}_{\mathrm{1}}, the momentum equation becomes:

(19){{\displaystyle \rho v}}_{,t}-{{\displaystyle \Sigma }}_{,x}+\beta {{\displaystyle \left({\displaystyle \frac{\sqrt{\mathrm{1}+\mathrm{2}{{\displaystyle \chi }}_{\mathrm{1}}\Theta }-\mathrm{1}}{{{\displaystyle \chi }}_{\mathrm{1}}}}\right)}}_{,x}=\rho f.

Let the vector of state variables (U) and fluxes (F) are:

U={\{v,\Sigma ,q,\Theta \}}^T,

(20)F={{\displaystyle \left\{{\displaystyle \frac{\mathrm{1}}{\rho }}\left(-\Sigma +\beta \left({\displaystyle \frac{\sqrt{\mathrm{1}+\mathrm{2}{{\displaystyle \chi }}_{\mathrm{1}}\Theta }-\mathrm{1}}{{{\displaystyle \chi }}_{\mathrm{1}}}}\right)\right),-{{\displaystyle \rho v}}_e^{\mathrm{2}},\bar{\alpha }\theta ,{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \tau }}_{\mathrm{0}}+{{\displaystyle \eta }}_{\mathrm{2}}}}\left({{\displaystyle \eta }}_{\mathrm{1}}{{\displaystyle k}}_{\mathrm{0}}\Theta +{{\displaystyle \eta }}_{\mathrm{2}}{{\displaystyle k}}_{\mathrm{0}}^{*}\left[\Theta +{\displaystyle \frac{({{\displaystyle \chi }}_{\mathrm{2}}-{{\displaystyle \chi }}_{\mathrm{1}})(\mathrm{1}+{{\displaystyle \chi }}_{\mathrm{1}}\Theta -\sqrt{\mathrm{1}+\mathrm{2}{{\displaystyle \chi }}_{\mathrm{1}}\Theta })}{{{\displaystyle \chi }}_{\mathrm{1}}^{\mathrm{2}}}}\right]\right),{\displaystyle \frac{{{\displaystyle v\rho v}}_e^{\mathrm{2}}}{\beta }}v+{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \rho c}}_{\mathrm{0}}}}q\right\}}}^T,

then, by considering Eqs. (19), (17), (18), and (16), the vector of source is:

(21)S={{\displaystyle \left\{f,\mathrm{0},-{\displaystyle \frac{{{\displaystyle \eta }}_{\mathrm{1}}}{{{\displaystyle \tau }}_{\mathrm{0}}+{{\displaystyle \eta }}_{\mathrm{2}}}}q,{\displaystyle \frac{\mathrm{1}}{{{\displaystyle c}}_{p\mathrm{0}}}}r\right\}}}^T,

and the Jacobian matrix becomes:

(22)A=\left(\begin{array}{llll}\mathrm{0}& -{\displaystyle \frac{\mathrm{1}}{\rho }}& \mathrm{0}& {\displaystyle \frac{\mathrm{2}\beta }{\rho \sqrt{\mathrm{1}+\mathrm{2}\Theta {{\displaystyle \chi }}_{\mathrm{1}}}}}\\ -{{\displaystyle \rho v}}_e^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {{\displaystyle \bar{\alpha }}}_{\mathrm{1}}\\ {\displaystyle \frac{{{\displaystyle v\rho v}}_e^{\mathrm{2}}}{\beta }}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{{{\displaystyle \rho c}}_{\mathrm{0}}}}& \mathrm{0}\end{array}\right),

where {{\displaystyle \bar{\alpha }}}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{{{\displaystyle \tau }}_{\mathrm{0}}+{{\displaystyle \eta }}_{\mathrm{2}}}}\left({{\displaystyle \eta }}_{\mathrm{1}}{{\displaystyle k}}_{\mathrm{0}}+{{\displaystyle \eta }}_{\mathrm{2}}{{\displaystyle k}}_{\mathrm{0}}^{*}\left[\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{{{\displaystyle \chi }}_{\mathrm{1}}}}({{\displaystyle \chi }}_{\mathrm{2}}-{{\displaystyle \chi }}_{\mathrm{1}})\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}}{\sqrt{\mathrm{1}+\mathrm{2}{{\displaystyle \chi }}_{\mathrm{1}}\Theta }}}\right)\right]\right). Eigenvalues of this Jacobian matrix are:

(23){{\displaystyle \bar{\lambda }}}_{\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}}=\pm {\displaystyle \frac{{{\displaystyle v}}_e}{\sqrt{\mathrm{2}}}}\sqrt{{{\displaystyle \bar{\alpha }}}_{\mathrm{2}}\pm \sqrt{{{\displaystyle \bar{\alpha }}}_{\mathrm{2}}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{4}{{\displaystyle \bar{\alpha }}}_{\mathrm{1}}{{\displaystyle v}}_t^{\mathrm{2}}}{k^{*}}}}},

where {{\displaystyle \bar{\alpha }}}_{\mathrm{2}}:={\displaystyle \frac{{{\displaystyle \bar{\alpha }}}_{\mathrm{1}}{{\displaystyle v}}_t^{\mathrm{2}}}{{{\displaystyle k}}_{\mathrm{0}}^{*}}}+{\displaystyle \frac{\mathrm{2}v}{\sqrt{\mathrm{2}{{\displaystyle \chi }}_{\mathrm{1}}\Theta +\mathrm{1}}}}+\mathrm{1}.

3 Multiresolution analysis (MRA) and grid adaptations

In this section, the concept of the MRA is presented by the wavelet theory. Then a family of interpolating wavelets is reviewed which has a simple computational algorithm with a physical meaning. Based on this wavelet transform, the concept and algorithm of a 1-D grid adaptation are presented.

3.1 MRA

Wavelets can simultaneously capture different local features of different supports (resolutions) of data. This capability makes them as mathematical microscopes. For f(x)\in L^{\mathrm{2}}(ℝ), wavelet transforms can decompose f(x) to a finite set of small and localized waves with different spatio–frequency supports. The wavelets, themselves, have two types: 1) functions \varphi (x) with low frequency content to estimate overall response; 2) functions \psi (x) with high frequency content to detect and represent localized fluctuations. The functions \varphi (x) and \psi (x) are known as the scale and wavelet functions, respectively. To cover all spatio-frequency components of data, dilated and shifted versions of \varphi (x) and \psi (x) are utilized in wavelet transforms. By dilation, frequency contents of the functions are adjusted and by shifting, different spatial locations can be studied by the functions. Intensity of variations can be measured by coefficients of the wavelet functions (\{\psi (x)\}): coefficient values are in accordance with variation magnitudes. This means, coefficients of large values concentrate automatically around considerable variations. The wavelet coefficients can then be considered as a criterion for detection of high gradient zones and grid adaptations.

3.2 Dubuc-Deslauriers (D-D) interpolating wavelets and wavelet transform

In this work, the D-D interpolating wavelets of order \mathrm{2}M-\mathrm{1} are used for grid adaptation. They have compact supports, as: Supp(\varphi )=[-\mathrm{2}M+\mathrm{1},\mathrm{2}M-\mathrm{1}], where \mathrm{2}M-\mathrm{1} is order of the interpolating scale function. In this family, both scaling and wavelet functions have the interpolation feature; their wavelets are defined as: \psi (x)=\varphi (\mathrm{2}x-\mathrm{1}) [61]. This wavelet family leads to a simple algorithm totally done in the physical domain; hence, transform coefficients have physical meaning. We survey the transform algorithm in the following.

Let us assume a dyadic grid as [62]:

(24){{\displaystyle V}}_j=\{{{\displaystyle x}}_{j,k}\in [\mathrm{0},\mathrm{1}]:{{\displaystyle x}}_{j,k}=k/{\mathrm{2}}^j\},j\in \mathbb{Z},k\in \{\mathrm{0},\mathrm{1},\mathrm{...},{\mathrm{2}}^j\},

where, j and k denote resolution levels (corresponding to resolution \mathrm{1}/{\mathrm{2}}^j) and spatial positions, respectively. Since {\displaystyle \frac{\mathrm{2}k}{{\mathrm{2}}^{j+\mathrm{1}}}}={\displaystyle \frac{k}{{\mathrm{2}}^j}} then {{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}k}={{\displaystyle x}}_{j,k}; this means {{\displaystyle V}}_j\subset {{\displaystyle V}}_{j+\mathrm{1}}: the core of MRA. This is known as the two-scale relationship in MRA. By repeating this procedure, it is easy to show that: {{\displaystyle V}}_{j-q}\subset \mathrm{...}{{\displaystyle V}}_{j-\mathrm{1}}\subset {{\displaystyle V}}_j\subset \mathrm{...}{{\displaystyle V}}_{j+p} for p,q\in \mathbb{N}. Let us assume two successive subspaces {{\displaystyle V}}_j and {{\displaystyle V}}_{j+\mathrm{1}}; points {{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}k+\mathrm{1}} belong to {{\displaystyle V}}_{j+\mathrm{1}} but are not in {{\displaystyle V}}_j. These points have the sampling step \mathrm{1}/{\mathrm{2}}^j and so their resolution level is j. Hence information on {{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}k+\mathrm{1}} can act as extra details to refine information from {{\displaystyle V}}_j into {{\displaystyle V}}_{j+\mathrm{1}}. So, the detail space {{\displaystyle W}}_j can be defined as:

(25){{\displaystyle W}}_j=\{{{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}k+\mathrm{1}}\in (\mathrm{0},\mathrm{1}):{{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}k+\mathrm{1}}=(\mathrm{2}k+\mathrm{1})/{\mathrm{2}}^{j+\mathrm{1}}\},j\in \mathbb{Z},k\in \{\mathrm{0},\mathrm{1},\mathrm{...},{\mathrm{2}}^j-\mathrm{1}\},

where {{\displaystyle V}}_j\oplus {{\displaystyle W}}_j={{\displaystyle V}}_{j+\mathrm{1}}.

Let us assume a continuous function f(x)\in L^{\mathrm{2}}(ℝ) defined on {{\displaystyle V}}_{J\max } (i.e., x\in {{\displaystyle V}}_{j\max }), where {{\displaystyle J}}_{\max } denotes the finest resolution level with the sampling step \mathrm{1}/{\mathrm{2}}^{J\max }. Considering the MRA, f(x) can be expanded in wavelet spaces, as [61]:

(26)f(x)={{\displaystyle \Sigma }}_{i=\mathrm{0}}^{{\mathrm{2}}^{J_{\min }}}{{\displaystyle c}}_{J_{\min },l}{{\displaystyle \varphi }}_{J_{\min },l}(x)+{{\displaystyle \Sigma }}_{j=J_{\min }}^{J_{\max }-\mathrm{1}}{{\displaystyle \Sigma }}_{n=\mathrm{0}}^{{\mathrm{2}}^j-\mathrm{1}}{{\displaystyle d}}_{j,n}{{\displaystyle \psi }}_{j,n}(x)={{\displaystyle Pf}}_{J_{\min }}+{{\displaystyle \Sigma }}_{j=J_{\min }}^{J_{\max }-\mathrm{1}}{{\displaystyle Qf}}_j,

where:{{\displaystyle J}}_{\min } is the coarsest resolution level; {{\displaystyle \varphi }}_{j,l}(x) and {{\displaystyle \psi }}_{j,n}(x) are dilated and shifted versions of \varphi (x) and \psi (x), respectively; i.e.: {{\displaystyle \varphi }}_{j,l}(x)=\varphi ({\mathrm{2}}^{j}x-l) & {{\displaystyle \psi }}_{j,n}(x)=\psi ({\mathrm{2}}^{j}x-n); {{\displaystyle c}}_{j,l} and {{\displaystyle d}}_{j,n} denote approximation and detail coefficients of the transform, respectively. Functions {{\displaystyle Pf}}_j:={{\displaystyle \Sigma }}_{n=\mathrm{0}}^{{\mathrm{2}}^j}{{\displaystyle c}}_{j,l}{{\displaystyle \varphi }}_{j,l}(x) and {{\displaystyle Qf}}_j:={{\displaystyle \Sigma }}_{n=\mathrm{0}}^{{\mathrm{2}}^j-\mathrm{1}}{{\displaystyle d}}_{j,n}{{\displaystyle \psi }}_{j,n}(x) show respectively the overall approximation and detail information of f(x) at resolution level j. To approximate f(x) on the successive finer resolution ({{\displaystyle Pf}}_{j+\mathrm{1}}), {{\displaystyle Qf}}_j should be added as extra information to {{\displaystyle Pf}}_j; this procedure can be repeated to increase resolution of approximated data. It should be mentioned that the dilated-shifted functions {{\displaystyle \varphi }}_{j,l}(x) and {{\displaystyle \psi }}_{j,n}(x) both have the interpolation property, as well.

In the D-D wavelets, coefficients of {{\displaystyle c}}_{j,l} and {{\displaystyle d}}_{j,n} can be evaluated regarding the interpolating feature. Due to this property (i.e.: \varphi (k)={{\displaystyle \delta }}_{k\mathrm{0}}\forall k\in \mathbb{Z}), the approximation coefficients ({{\displaystyle c}}_{J_{\min },l}) are equal to {{\displaystyle c}}_{J_{\min },l}=f({{\displaystyle x}}_{J_{\min },l}): sampled values of f(x) at points {{\displaystyle x}}_{J_{\min },l}\in {{\displaystyle V}}_{J_{\min }}. And this means that {{\displaystyle Pf}}_j({{\displaystyle x}}_{J_{\min },l})=f({{\displaystyle x}}_{J_{\min },l}).

Let assume {{\displaystyle Pf}}_j({{\displaystyle x}}_{j,l})=f({{\displaystyle x}}_{j,l}); regarding the two-scale relationship, it is clear: {{\displaystyle x}}_{j,l}={{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l} and so {{\displaystyle Pf}}_j({{\displaystyle x}}_{j,l})={{\displaystyle Pf}}_j({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l})=f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l}). This relationship, however, do not valid for odd-numbered grid points {{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}}\in {{\displaystyle V}}_{j+\mathrm{1}}; that is: {{\displaystyle Pf}}_j({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}})\ne f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}}). The difference between f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}}) and their projections at {{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}}, {{\displaystyle Pf}}_j({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}}), are defined as the detail coefficients [70]. Having these differences and {{\displaystyle Pf}}_j({{\displaystyle x}}_{j,l}) values, the successive finer approximation {{\displaystyle Pf}}_{j+\mathrm{1}} can be retrieved, as: {{\displaystyle Pf}}_{j+\mathrm{1}}={{\displaystyle Pf}}_j+{{\displaystyle Qf}}_j [62].

In the D-D wavelets, {{\displaystyle Pf}}_j({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}}) can be obtained by a local Lagrange interpolation. For this interpolation, for the D-D wavelet of order \mathrm{2}M-\mathrm{1}, \mathrm{2}M most neighbor points \{{{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{2}n}\} are selected around {{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}}, where n\in \{-M+\mathrm{1},-M+\mathrm{2},\mathrm{...},M\}. It is easy to show that for case M=\mathrm{2}, {{\displaystyle Pf}}_j({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}}) can be obtained as:

(27){{\displaystyle Pf}}_j({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}})=\left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{16}}}\right)f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l-\mathrm{2}})+\left({\displaystyle \frac{\mathrm{9}}{\mathrm{16}}}\right)f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l})+\left({\displaystyle \frac{\mathrm{9}}{\mathrm{16}}}\right)f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{2}})+\left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{16}}}\right)f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{4}})

And then the detail coefficients are:

(28){{\displaystyle d}}_{j,l}=f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}})-{{\displaystyle Pf}}_j({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}l+\mathrm{1}}).

For grid points of finite length, the boundary wavelets, introduced by Donoho [71] are used in the vicinity of edge points. For grid points \{\{{{\displaystyle x}}_i\};i=\{\mathrm{0},\mathrm{1},\mathrm{...},{\mathrm{2}}^{J_{\max }}\}\} and for case M=\mathrm{2}, it is easy to show that {{\displaystyle Pf}}_j({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}n+\mathrm{1}}) are:

(29){{\displaystyle Pf}}_j({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}n+\mathrm{1}})={\displaystyle \frac{\mathrm{1}}{\mathrm{16}}}\{\begin{array}{l}\mathrm{5}f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{0}})+\mathrm{15}f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}})-\mathrm{5}f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{4}})+f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{6}}),\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}n=\mathrm{0},\\ -f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}n-\mathrm{2}})+\mathrm{9}f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}n})+\mathrm{9}f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}n+\mathrm{2}})-f({{\displaystyle x}}_{j+\mathrm{1},\mathrm{2}n+\mathrm{4}}),\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}n=\mathrm{1},\mathrm{...},{{\displaystyle N}}_j-\mathrm{2}\\ f({{\displaystyle x}}_{j+\mathrm{1},{{\displaystyle N}}_j-\mathrm{6}})-\mathrm{5}f({{\displaystyle x}}_{j+\mathrm{1},{{\displaystyle N}}_j-\mathrm{4}})+\mathrm{15}f({{\displaystyle x}}_{j+\mathrm{1},{{\displaystyle N}}_j-\mathrm{2}})+\mathrm{5}f({{\displaystyle x}}_{j+\mathrm{1},{{\displaystyle N}}_j}),\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}n={{\displaystyle N}}_j-\mathrm{1}\end{array},

where {{\displaystyle N}}_j={\mathrm{2}}^j and j=\{{{\displaystyle J}}_{\min },{{\displaystyle J}}_{\min }+\mathrm{1},\mathrm{...},{{\displaystyle J}}_{\max }-\mathrm{1}\}.

3.3 One dimensional grid adaptations by MRA

As mentioned before, the adaptation criterion is magnitude of wavelet (detail) coefficients. By defining a predefined thresholds {{\displaystyle \in }}_j ({{\displaystyle \in }}_j can be different for each resolution level j), detail coefficients can be updated as:

(30){{\displaystyle \widehat{d}}}_{j,n}=\{\begin{array}{l}\mathrm{0},if\left|{{\displaystyle d}}_{j,n}\right|\langle {{\displaystyle \in }}_j,\\ {{\displaystyle d}}_{j,n},if\left|{{\displaystyle d}}_{j,n}\right|\ge {{\displaystyle \in }}_j.\end{array}

By this thresholding, wavelet coefficient and corresponding grid points are omitted. Holmström [69] showed that such truncation error is bounded and approaches to zero as {{\displaystyle \in }}_j\to \mathrm{0}.

This thresholding can be done independently for each resolution j where j=\{{{\displaystyle J}}_{\min },\mathrm{...},{{\displaystyle J}}_{\max }-\mathrm{1}\} with either constant ({{\displaystyle \in }}_j=\in) or level-dependent ({{\displaystyle \in }}_j) thresholds. For the constant-value thresholding, using of normalized detail coefficients is recommended, as: {{\displaystyle d}}_{j,n}={{\displaystyle d}}_{j,n}/{{\displaystyle f}}_j^{\mathrm{r}\mathrm{e}\mathrm{t}}. Here, the first approach with the normalized factor {{\displaystyle f}}_j^{\mathrm{r}\mathrm{e}\mathrm{f}}={{\displaystyle \max }}_{j}f({{\displaystyle x}}_j) is used. In practical problems, the recommended value of {{\displaystyle \in }}_j is around {\mathrm{10}}^{-\mathrm{4}}; suitable value can be adjusted by the trial-error method.

Aforementioned adaptation procedure is for scalar functions. For a system of vector form, the resultant adapted grid is superposition of all adapted grids where each one is obtained independently. By this, effects of all solutions are considered in the adaptation procedure.

3.4 MRA-based modification of adapted grids

To prevent abrupt variations in adapted grids and so to preserve the numerical stability, a rechecking of adapted grids is essential. Here, this stage is also performed by the MRA: for each adapted grid at the resolution level j, it is essential to control existence of sufficient surrounding grid points at two successive resolution levels: i) in the same resolution level j: to have sufficient number of surrounding grids: i.e., {{\displaystyle N}}_s points at each side, ii) the successive coarser scale j-\mathrm{1}: there are enough number of adapted points for resolution level j-\mathrm{1}; i.e., {{\displaystyle N}}_c points at each side in the successive resolution level j-\mathrm{1}. These controlling/adding procedures are started from the finest resolution {{\displaystyle J}}_{\max }-\mathrm{1} and repeated until reaching the level {{\displaystyle J}}_{\min }+\mathrm{1}. At the coarsest level {{\displaystyle J}}_{\min }, only the first checking/adding step (in the same resolution) is done (see also Appendix A).

4 Reconstructions by average interpolating radial basis functions (RBFs)

4.1 Average interpolating RBFs

For a cell {{\displaystyle I}}_j and corresponding symmetric neighbor cells, a stencil {{\displaystyle S}}_j is defined as:

{{\displaystyle S}}_j=\{{{\displaystyle I}}_{j-q},\mathrm{...},{{\displaystyle I}}_{j-\mathrm{1}},{{\displaystyle I}}_{j+\mathrm{1}},\mathrm{...},{{\displaystyle I}}_{j+q}\},where size of the stencil {{\displaystyle S}}_j is Q=\mathrm{2}q+\mathrm{1}.

An average interpolating RBF on the set {{\displaystyle S}}_j can be defined as [58]:

(32)s^j(x):={{\displaystyle \Sigma }}_{k=j-q}^{j+q}{{\displaystyle c}}_k[{{\displaystyle \mu }}_k^y\varphi (\Vert x-y\Vert )]+{{\displaystyle p}}_n(x),x\in {{\displaystyle S}}_j,

where: \varphi (x) denotes a RBF; {{\displaystyle c}}_k are unknown coefficients; {{\displaystyle p}}_n(x) shows a polynomial of order at most n (or of degree at most n-\mathrm{1}); \Vert \cdot \Vert denotes the Euclidean norm; and {\mu }_i^{y}Z:=\frac{1}{\left|I_i\right|}{\int }_{I_i}Z\left(y\right)dy is average of Z with respect to the variable y over the cell {{\displaystyle I}}_i, in which \mu is called the averaging operator and for one-dimensional elements, \left|{{\displaystyle I}}_i\right|:={{\displaystyle x}}_{i+\mathrm{1}/\mathrm{2}}-{{\displaystyle x}}_{i-\mathrm{1}/\mathrm{2}} is length of the cell {{\displaystyle I}}_i ({{\displaystyle x}}_{i\pm \mathrm{1}/\mathrm{2}} are cell edges). Regarding this definition of the averaging operator, {{\displaystyle \mu }}_k^y\varphi (\Vert x-y\Vert ) yields:

(33){\mu }_k^y\phi (\Vert x-y\Vert ):=\frac{1}{\left|I_k\right|}{\int }_{I_k}\phi (\Vert x-y\Vert )dy,

The polynomial order (n) is determined by the order of the RBF \varphi. Some RBFs and corresponding parameters are presented in Table 1; there, d denotes the dimension of RBFs. where c and \alpha are the shape parameters, controlling accuracy and stability of RBF-based interpolations. In this study, \alpha is assumed to be \mathrm{1}/\mathrm{2} for the MQ-RBFs and -\mathrm{1}/\mathrm{2} for the IMQ-RBFs. Different choosing approaches of c are reviewed in subsection 4.3.

In Eq. (32), Q+n unknown parameters should be determined; Q relationships can be obtained by the reconstruction requirements:

(34){{\displaystyle \mu }}_k^{x}s={{\displaystyle \bar{u}}}_k,

where {{\displaystyle \bar{u}}}_k denotes the average of solution u(x):=u(x,t) on the cell {{\displaystyle I}}_k. Remaining unknown coefficients can be determined by the constraints [58]:

(35){{\displaystyle \Sigma }}_k{{\displaystyle c}}_k{{\displaystyle \mu }}_k^x(x^{\alpha })=\mathrm{0},\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}all\hspace{0.17em}\mathrm{0}\le \alpha \langle n.

These equations and constraints yield a (Q+n)\times (Q+n) linear system:

(36)\left[\begin{array}{ll}\bar{A}& \bar{P}\\ {{\displaystyle \bar{P}}}^T& \mathrm{0}\end{array}\right]\left[\begin{array}{l}c\\ d\end{array}\right]=\left[\begin{array}{l}\bar{u}\\ \mathrm{0}\end{array}\right],

where: {{\displaystyle \overline{A}}}_{i,l}={{\displaystyle \mu }}_i^x{{\displaystyle \mu }}_l^y\phi (\Vert x-y\Vert ) for i,l=\{j-q,\mathrm{...},j+q\}; {{\displaystyle \bar{P}}}_{i,l}={{\displaystyle \mu }}_i^x(x^{l-\mathrm{1}}) for i=\{j-q,\mathrm{...},j+q\} & i=\{\mathrm{1},\mathrm{2},\mathrm{...},n\}; {{\displaystyle c}}_i={{\displaystyle c}}_i for i=\{j-q,\mathrm{...},j+q\}; {{\displaystyle d}}_i={{\displaystyle d}}_i for i=\{\mathrm{1},\mathrm{2},\mathrm{...}n\} in which {{\displaystyle d}}_i is coefficient of the term x^{i-\mathrm{1}} (in the {{\displaystyle p}}_n(x)); and {{\displaystyle \bar{u}}}_i={{\displaystyle \bar{u}}}_i denotes the average solution on the cell {{\displaystyle I}}_i for i=\{j-q,\mathrm{...},j+q\}.

Having functions s^j(x) defined on sets {{\displaystyle S}}_j, a piecewise function can be defined as:

(37)s(x):={{\displaystyle \mathrm{\Sigma }}}_k{s}^k(x){{\displaystyle \tilde{s}}}_k(x),

where {{\displaystyle \tilde{s}}}_k(x) denotes the unit function in which {{\displaystyle \tilde{s}}}_k(x)=\mathrm{1} for x\in {{\displaystyle I}}_j and is zero elsewhere.

4.2 Monotonic piecewise interpolations with RBFs

To guarantee monotonicity across cell edges, reconstructions are rescaled as [38]:

(38){{\displaystyle s}}_m^j(x)={{\displaystyle \bar{u}}}_j+{{\displaystyle \theta }}_j\left(s^j(x)-{{\displaystyle \bar{u}}}_j\right),\mathrm{0}\le {{\displaystyle \theta }}_j\le \mathrm{1},

where {{\displaystyle \theta }}_j is a nonlinear scaling limiter, defined as:

(39){{\displaystyle \theta }}_j:=\{\begin{array}{l}\min \left({\displaystyle \frac{{{\displaystyle M}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}-{{\displaystyle \bar{u}}}_j}{{{\displaystyle M}}_j-{{\displaystyle \bar{u}}}_j}}+{\displaystyle \frac{{{\displaystyle M}}_{j-\frac{\mathrm{1}}{\mathrm{2}}}-{{\displaystyle \bar{u}}}_j}{{{\displaystyle M}}_j-{{\displaystyle \bar{u}}}_j}},\mathrm{1}\right),if{{\displaystyle \bar{u}}}_{j-\mathrm{1}}\langle {{\displaystyle \bar{u}}}_j\langle {{\displaystyle \bar{u}}}_{j+\mathrm{1}},\\ \min \left({\displaystyle \frac{{{\displaystyle M}}_{j-\frac{\mathrm{1}}{\mathrm{2}}}-{{\displaystyle \bar{u}}}_j}{{{\displaystyle M}}_j-{{\displaystyle \bar{u}}}_j}}+{\displaystyle \frac{{{\displaystyle M}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}-{{\displaystyle \bar{u}}}_j}{{{\displaystyle M}}_j-{{\displaystyle \bar{u}}}_j}},\mathrm{1}\right)if{{\displaystyle \bar{u}}}_{j-\mathrm{1}}\rangle {{\displaystyle \bar{u}}}_j\rangle {{\displaystyle \bar{u}}}_{j+\mathrm{1}},\\ \mathrm{1},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}otherwise,\end{array}

in which:

{{\displaystyle M}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}=\max \left({\displaystyle \frac{{{\displaystyle \bar{u}}}_{j+\mathrm{1}}+{{\displaystyle \bar{u}}}_j}{\mathrm{2}}},s^{j+\mathrm{1}}\left({{\displaystyle x}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}\right)\right),{{\displaystyle M}}_{j-\frac{\mathrm{1}}{\mathrm{2}}}=\max \left({\displaystyle \frac{{{\displaystyle \bar{u}}}_j+{{\displaystyle \bar{u}}}_{j-\mathrm{1}}}{\mathrm{2}}},s^{j-\mathrm{1}}\left({{\displaystyle x}}_{j-\frac{\mathrm{1}}{\mathrm{2}}}\right)\right),

(40){{\displaystyle m}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}=\min \left({\displaystyle \frac{{{\displaystyle \bar{u}}}_{j+\mathrm{1}}+{{\displaystyle \bar{u}}}_j}{\mathrm{2}}},s^{j+\mathrm{1}}\left({{\displaystyle x}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}\right)\right),{{\displaystyle m}}_{j-\frac{\mathrm{1}}{\mathrm{2}}}=\min \left({\displaystyle \frac{{{\displaystyle \bar{u}}}_j+{{\displaystyle \bar{u}}}_{j-\mathrm{1}}}{\mathrm{2}}},s^{j-\mathrm{1}}\left({{\displaystyle x}}_{j-\frac{\mathrm{1}}{\mathrm{2}}}\right)\right),

and:

(41){{\displaystyle M}}_j=\underset{x\in {{\displaystyle I}}_j}{{\max s}^j}(x),{{\displaystyle m}}_j=\underset{x\in {{\displaystyle I}}_j}{{\min s}^j}(x).

For monotonic solutions, parameters {{\displaystyle M}}_j and {{\displaystyle m}}_j will be:

(42){{\displaystyle M}}_j=\max \left(s^j\left({{\displaystyle x}}_{j-\frac{\mathrm{1}}{\mathrm{2}}}\right),s^j\left({{\displaystyle x}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}\right)\right),{{\displaystyle m}}_j=\min \left(s^j\left({{\displaystyle x}}_{j-\frac{\mathrm{1}}{\mathrm{2}}}\right),s^j\left({{\displaystyle x}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}\right)\right).

This rescaling is done in a way that reconstructions remain conservative, i.e., {{\displaystyle \mu }}_j^x{{\displaystyle s}}_m^j(x)={{\displaystyle \mu }}_j^x{s}^j(x). Using this modified function {{\displaystyle s}}_m^j(x), the monotonized piecewise reconstruction becomes:

(43)s(x):={{\displaystyle \mathrm{\Sigma }}}_k{{\displaystyle s}}_m^k(x){{\displaystyle \tilde{s}}}_k(x)

4.3 On choosing a good shape parameter c

It is known that there is a trade-off between good approximation feature and good numerical stability for interpolations obtained by RBFs [60]. Obtaining theoretically achievable accuracy is impossible due to numerically instabilities. This drawback is because of a large condition number of the interpolation matrix in Eq. (36). Such trade-off relationship is referred to the uncertainty relation or the trade-off principle. It is also shown that it is impossible to keep small both the numerical error and the condition number, simultaneously [60]. The shape parameter c has considerable effects on RBFs performance. It controls trade-off between interpolation accuracy and instability.

For handling the trade-off principle for finding an optimum value for the shape parameter, several ad-hocsolutions were proposed; most of them were based on some numerical sensitive analysis by different test problems [52,72‒74]; for a general review, see also [75,76]. Several algorithms have been proposed for treatment of the uncertainty relation; some of which are:

1) Utilizing of indirect solvers for achieving stable results [77‒80]: this overcomes drawbacks of direct solvers (with large condition numbers). The main shortcoming of these approaches may be their computation cost,

2) Integration of the regularization concept with RBF-based solvers to attain stable solutions [81,82],

3) Following the preconditioning technique to decrease condition numbers [83],

4) Employing of spatially variable shape parameters [84,85],

5) Using of an optimization approach [86‒89].

It should be mentioned that above mentioned approaches may also be integrated with each other.

In this work, we follow the simple strategy explained in [52,74]. The main idea is to perform a sensitivity analysis by a series of RBF-based interpolations to determine an optimum or a good value of the shape parameter. This study is done for different cell-sizes ({{\displaystyle \Delta x}}_i) and different values of the shape parameters, simultaneously. In general, the optimal value of the shape parameter can depend on RBFs and also features of functions f(x) needed to be interpolated (i.e. smoothness or discontinuity of f(x)).

For simplicity, we only consider 1-D functions; for smooth-high gradient functions, f(x) is recommended to be as follows [52]:

(44){{\displaystyle f}}_{\mathrm{1}}:={{\displaystyle f}}_{\mathrm{1}}(x)=\sin (\mathrm{2}\pi x),x\in [\mathrm{0},\mathrm{1}].

In this study, also, for solutions with discontinuities, performance of RBFs is controlled for the unit step function, defined as:

(45){{\displaystyle f}}_{\mathrm{2}}:={{\displaystyle f}}_{\mathrm{2}}(x)=H(x-\mathrm{0.5}),x\in [\mathrm{0},\mathrm{1}],

where H(x) is the Heaviside function.

For the sensitivity analysis, the domain \Omega =[\mathrm{0},\mathrm{1}] is discretized uniformly with different sampling steps {{\displaystyle \Delta x}}_i=\mathrm{1}/{\mathrm{2}}^i for i=\mathrm{5},\mathrm{6},\mathrm{7},\mathrm{8},\mathrm{9}. At cell-edges {{\displaystyle x}}_{j+\mathrm{1}/\mathrm{2}}, reconstructed values {{\displaystyle u}}_{j+\mathrm{1}/\mathrm{2}}^L:={{\displaystyle s}}_m^j({{\displaystyle x}}_{j+\mathrm{1}/\mathrm{2}}) and {{\displaystyle u}}_{j+\mathrm{1}/\mathrm{2}}^R:={{\displaystyle s}}_m^{j+\mathrm{1}}({{\displaystyle x}}_{j+\mathrm{1}/\mathrm{2}}) should be evaluated by the monotonized piecewise interpolations (Eq. (43)); for this reason, proper error estimators could be those considering reconstruction errors at cell-interfaces. So, the error definition proposed in [52,90] is utilized, as:

(46)Err:={{\displaystyle \Vert {{\displaystyle u}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{-}-f\left({{\displaystyle x}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}\right)\Vert }}_{{{\displaystyle L}}_{\Delta x}^p}+{{\displaystyle \Vert {{\displaystyle u}}_{i+\mathrm{1}/\mathrm{2}}^{+}-f\left({{\displaystyle x}}_{i+\frac{\mathrm{1}}{\mathrm{2}}}\right)\Vert }}_{{{\displaystyle L}}_{\Delta x}^p},

where \Delta x=\mathrm{1}/N and the norm is:

(47){{\displaystyle \Vert {{\displaystyle (\cdot )}}_i\Vert }}_{{{\displaystyle L}}_{\Delta x}^p}:={({{\displaystyle \Sigma }}_{i=\mathrm{1}}^{N-\mathrm{1}}{\left|{{\displaystyle (\cdot )}}_i\right|}^p)}^{\mathrm{1}/p}for\hspace{0.17em}p\langle \infty ,and\hspace{0.17em}{{\displaystyle \Vert {{\displaystyle (\cdot )}}_i\Vert }}_{{{\displaystyle L}}_{\Delta x}^{\infty }}:={{\displaystyle \max }}_i\left|{{\displaystyle (\cdot )}}_i\right|.

Results are presented for the functions {{\displaystyle f}}_{\mathrm{1}} and {{\displaystyle f}}_{\mathrm{2}} in Fig. 1; there G-, MQ-, and IMQ-RBFs are considered where p=\mathrm{1}. The results offer that:

1) In general, all of the RBFs are sensitive for both small c and small \Delta x values,

2) The optimum value of c can depend on \Delta x,

3) The G-RBF is not so sensitive for \Delta x and c (in comparison with MQ and IMQ RBFs),

4) For the smooth function {{\displaystyle f}}_{\mathrm{1}}, for \Delta x\ge \mathrm{1}/{\mathrm{2}}^{\mathrm{8}}, optimum values of the shape parameters c are: i) for\hspace{0.17em}G-RBF: {{\displaystyle c}}_{opt}=\mathrm{6.5}; ii) for MQ-RBF: {{\displaystyle c}}_{opt}=\mathrm{13}; iii) for IMQ-RBF: {{\displaystyle c}}_{opt}=\mathrm{4.5},

5) For {{\displaystyle f}}_{\mathrm{1}}, however, these optimum values, could change for smaller sampling steps, i.e., for \Delta x=\mathrm{1}/{\mathrm{2}}^{\mathrm{9}}. For IMQ-RBF, this sensitivity is high: the optimum value {{\displaystyle c}}_{opt}=\mathrm{4.5} leads to large errors for the test function {{\displaystyle f}}_{\mathrm{2}},

6) For the test function {{\displaystyle f}}_{\mathrm{2}}, MQ and IMQ RBFs are very sensitive for small values of \Delta x and c. To decrease this sensitivity, large-enough values of c with smaller condition numbers are recommended; see results for \Delta x=\mathrm{1}/{\mathrm{2}}^{\mathrm{9}}.

In conclusion, proper choosing of the shape parameter for the MQ- and IMQ-RBFs needs especial attention. For these RBFs, selected values of the shape parameters are: 1) for MQ-RBF: c={{\displaystyle c}}_{opt}=\mathrm{13} (with small-enough condition number); 2) for IMQ-RBF: c=\mathrm{6} (a larger value than {{\displaystyle c}}_{opt} to decrease the condition number). For the G-RBF, some simulations will be done with a small c value (c=\mathrm{1}), the optimum value ({{\displaystyle c}}_{opt}=\mathrm{6.5}) and a larger one (c=\mathrm{10}) to study effects of the shape parameter.

5 RBF-based higher order central high resolution schemes

In this section, it is shown how to integrate the monotonized piecewise RBF-based reconstructions with standard higher order central high resolution schemes. For this purpose, relationship between reconstructions obtained by polynomials and RBFs are studied. Based on this, the semi-discrete form of RBF-based central high resolution scheme is presented. High resolution schemes refer to methods using higher order approximation over smooth solutions and limited/lower-order ones around discontinuities. At the end of this section, the computational complexity of the RBF-based central high resolution scheme is investigated. Then, computation cost of this method is compared with other third-order central high resolution schemes.

5.1 Relationship between RBF- and polynomial-based reconstructions

The original high-order high resolution schemes use polynomials in their reconstruction stages; based on this type of reconstruction, corresponding fully-discrete and semi-discrete forms were provided. In this study, RBFs are used in the reconstruction stage. Hence, relationship between RBFs and polynomials is important for proper using of these fully or semi-discrete forms. In the following such relationships are studied for the G, MQ and IMQ RBFs; in all cases it is assumed {{\displaystyle p}}_n(x)=\mathrm{0}.

5.1.1 G-RBFs

It is easy to show that the averaging operator {{\displaystyle \mu }}_i^y\phi (\Vert x-y\Vert ) on the cell {{\displaystyle I}}_i is:

(48){{\displaystyle \mu }}_i^y\phi (\Vert x-y\Vert )={\displaystyle \frac{-\sqrt{\pi }}{\mathrm{2}\sqrt{c}{{\displaystyle \Delta x}}_i}}[Erf(\sqrt{c}(x-{{\displaystyle x}}_{i+\mathrm{1}/\mathrm{2}}))-Erf],

where Erf(x):={\displaystyle \frac{\mathrm{2}}{\sqrt{\pi }}}{{\displaystyle \int }}_{\mathrm{0}}^x{e}^{{-t}^{\mathrm{2}}}dt denotes the error function, {{\displaystyle \Delta x}}_i={{\displaystyle x}}_{i+\mathrm{1}/\mathrm{2}}-{{\displaystyle x}}_{i-\mathrm{1}/\mathrm{2}} is the cell length and {{\displaystyle x}}_i is the cell center. For uniform cells, using the Taylor series about {{\displaystyle x}}_i leads:

(49){{\displaystyle \mu }}_i^y\varphi (\Vert x-y\Vert )={\displaystyle \frac{\sqrt{\pi }Erf\left[{\displaystyle \frac{\sqrt{c}{{\displaystyle \Delta x}}_i}{\mathrm{2}}}\right]}{\sqrt{c}{{\displaystyle \Delta x}}_i}}-{ce}^{-\frac{{c{{\displaystyle \Delta x}}_i}^{\mathrm{2}}}{\mathrm{4}}}{(x-{{\displaystyle x}}_i)}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}\left({c^{\mathrm{2}}e}^{-\frac{{c{{\displaystyle \Delta x}}_i}^{\mathrm{2}}}{\mathrm{4}}}(-\mathrm{6}+{c{{\displaystyle \Delta x}}_i}^{\mathrm{2}})\right){(x-{{\displaystyle x}}_i)}^{\mathrm{4}}+O({{{\displaystyle \Delta x}}_i}^{\mathrm{6}}).

5.1.2 MQ-RBFs

In this case, for \alpha =\mathrm{1}/\mathrm{2}, the averaging operator {{\displaystyle \mu }}_i^y\phi (\Vert x-y\Vert ) on the cell {{\displaystyle I}}_i becomes:

(50){{\displaystyle \mu }}_i^y\phi (\Vert x-y\Vert )={\displaystyle \frac{-\mathrm{1}}{\mathrm{2}\sqrt{c}{{\displaystyle \Delta x}}_i}}(Z(x-{{\displaystyle x}}_{i-\mathrm{1}/\mathrm{2}})),

where Z(x):=x\sqrt{c(\mathrm{1}+{cx}^{\mathrm{2}})}+\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}[\sqrt{c}x] and {{\displaystyle x}}_{i\pm \mathrm{1}/\mathrm{2}} denotes cell edges. On uniform cells, for c=\mathrm{1}, expanding {{\displaystyle \mu }}_i^y\phi (\Vert x-y\Vert ) by the Taylor series around {{\displaystyle x}}_i yields:

(51){{\displaystyle \mu }}_i^y\varphi (\Vert x-y\Vert )=\left({\displaystyle \frac{\sqrt{\mathrm{4}+{{\displaystyle \Delta x}}_i^{\mathrm{2}}}}{\mathrm{4}}}+{\displaystyle \frac{arcsinh\left[{\displaystyle \frac{{{\displaystyle \Delta x}}_i}{\mathrm{2}}}\right]}{{{\displaystyle \Delta x}}_i}}\right)+{\displaystyle \frac{{(x-{{\displaystyle x}}_i)}^{\mathrm{2}}}{{(\mathrm{4}+{{\displaystyle \Delta x}}_i^{\mathrm{2}})}^{\mathrm{1}/\mathrm{2}}}}-{\displaystyle \frac{\mathrm{4}{(x-{{\displaystyle x}}_i)}^{\mathrm{4}}}{{(\mathrm{4}+{{\displaystyle \Delta x}}_i^{\mathrm{2}})}^{\mathrm{5}/\mathrm{2}}}}+O({{{\displaystyle \Delta x}}_i}^{\mathrm{6}}).

5.1.3 IMQ-RBFs

Let us assume \alpha =-\mathrm{1}/\mathrm{2}, then the average interpolating operator {{\displaystyle \mu }}_i^y\phi (\Vert x-y\Vert ) on the cell {{\displaystyle I}}_i is:

(52){{\displaystyle \mu }}_i^y\phi (\Vert x-y\Vert )={\displaystyle \frac{\mathrm{1}}{\sqrt{c}{{\displaystyle \Delta x}}_i}}(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}[\sqrt{c}(x-{{\displaystyle x}}_{i+\mathrm{1}/\mathrm{2}})]-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}[\sqrt{c}(x-{{\displaystyle x}}_{i-\mathrm{1}/\mathrm{2}})]),

using the Taylor series about {{\displaystyle x}}_i on uniform cells yields:

(53){{\displaystyle \mu }}_i^y\phi (\Vert x-y\Vert )={\displaystyle \frac{\mathrm{2}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\left[{\displaystyle \frac{\sqrt{c}{{\displaystyle \Delta x}}_i}{\mathrm{2}}}\right]}{\sqrt{c}{{\displaystyle \Delta x}}_i}}-{\displaystyle \frac{\mathrm{4}c{(x-{{\displaystyle x}}_i)}^{\mathrm{2}}}{{(\mathrm{4}+{{{\displaystyle c\Delta x}}_i}^{\mathrm{2}})}^{\mathrm{3}/\mathrm{2}}}}-{\displaystyle \frac{{\mathrm{8}c}^{\mathrm{2}}(-\mathrm{6}+{{{\displaystyle c\Delta x}}_i}^{\mathrm{2}}){(x-{{\displaystyle x}}_i)}^{\mathrm{4}}}{{(\mathrm{4}+{{{\displaystyle c\Delta x}}_i}^{\mathrm{2}})}^{\mathrm{7}/\mathrm{2}}}}+O({{{\displaystyle \Delta x}}_i}^{\mathrm{6}}).

Regarding Eqs. (49), (51) and (53), it is clear that the RBF-based reconstructions can approximate a parabolic polynomial with the truncation error of the order O({{{\displaystyle \Delta x}}_i}^{\mathrm{4}}). This reveals if these RBFs are used for the reconstruction stage, then corresponding central high resolution schemes have at least third order accuracy: high order central schemes (such accuracy order also was shown numerically for WENO-based RBF reconstructions with semi-uniform unstructured cells in [58]).

5.2 Semi-discrete form of central high resolution schemes

Using the Taylor series, it is shown that the RBF-based average interpolations have a close relationship with polynomial-based interpolation. In different RBFs, corresponding interpolation can be approximated with a parabolic polynomial with the error of order O({{{\displaystyle \Delta x}}_i}^{\mathrm{4}}). Based on this approximation, semi-discrete forms of polynomial-based central high resolution schemes can be used. The basic formulation is derived by the reconstruction-evolution-projection concept. For the balance equation {{\displaystyle u}}_{,t}+F{{\displaystyle (u)}}_{,x}=s(x,t) (where u:=u(x,t) and s:s(x,t) is the source term), the semi-discrete form is [38]:

(54){\displaystyle \frac{d{{\displaystyle \bar{u}}}_j(t)}{dt}}+{\displaystyle \frac{{{\displaystyle F}}_{j+\mathrm{1}/\mathrm{2}}^{*}-{{\displaystyle F}}_{j-\mathrm{1}/\mathrm{2}}^{*}}{{{\displaystyle \Delta x}}_j}}={{\displaystyle s}}_j,

where {{\displaystyle F}}_{j\pm \mathrm{1}/\mathrm{2}}^{*}=F^{*}({{\displaystyle u}}_{j\pm \mathrm{1}/\mathrm{2}}) are numerical fluxes defined as:

(55){{\displaystyle F}}_{j\pm \mathrm{1}/\mathrm{2}}^{*}:={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(F\left({{\displaystyle u}}_{j\pm \frac{\mathrm{1}}{\mathrm{2}}}^R(t)\right)+F\left({{\displaystyle u}}_{j\pm \frac{\mathrm{1}}{\mathrm{2}}}^L(t)\right)\right)-\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{{\displaystyle a}}_{j\pm \frac{\mathrm{1}}{\mathrm{2}}}\right)\left({{\displaystyle u}}_{j\pm \frac{\mathrm{1}}{\mathrm{2}}}^R(t)-{{\displaystyle u}}_{j\pm \frac{\mathrm{1}}{\mathrm{2}}}^L(t)\right),

where: u and f(u) denote the state variable and flux, respectively; {{\displaystyle \bar{u}}}_j is the average solution on the cell {{\displaystyle I}}_j; {{\displaystyle u}}_{j+\mathrm{1}/\mathrm{2}}^L and {{\displaystyle u}}_{j+\mathrm{1}/\mathrm{2}}^R show respectively the left and right reconstructed values of u at the cell-edge {{\displaystyle x}}_{j+\mathrm{1}/\mathrm{2}} (that is:

{{\displaystyle u}}_{j+\mathrm{1}/\mathrm{2}}^L:={{\displaystyle s}}_m^j\left({{\displaystyle x}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}\right)

and

{{\displaystyle u}}_{j+\mathrm{1}/\mathrm{2}}^R:={{\displaystyle s}}_m^{j+\mathrm{1}}\left({{\displaystyle x}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}\right)

); {{\displaystyle \Delta x}}_j:={{\displaystyle x}}_{j+\mathrm{1}/\mathrm{2}}-{{\displaystyle x}}_{j-\mathrm{1}/\mathrm{2}} denotes length of the cell {{\displaystyle I}}_j; and the parameter {{\displaystyle a}}_{J\pm \frac{\mathrm{1}}{\mathrm{2}}} denotes the maximum value of absolute values of eigenvalues at the cell-edges {{\displaystyle x}}_{j\pm \mathrm{1}/\mathrm{2}}. The eigenvalues are evaluated from the Jacobian matrix of the hyperbolic system.

For the time integration, the third-order TVD Runge-Kutta method is used [91]. The time-stepping from t^n to t^{n+\mathrm{1}} can be summarized as:

{{\displaystyle u}}_j^{(\mathrm{1})}={{\displaystyle u}}_j^{(\mathrm{0})}+\Delta tC({{\displaystyle u}}_j^{(\mathrm{0})}),

(56){{\displaystyle u}}_j^{(\mathrm{2})}={{\displaystyle u}}_j^{(\mathrm{1})}+{\displaystyle \frac{\Delta t}{\mathrm{4}}}[-\mathrm{3}C({{\displaystyle u}}_j^{(\mathrm{0})})+C({{\displaystyle u}}_j^{(\mathrm{1})})+\mathrm{8}C({{\displaystyle u}}_j^{(\mathrm{2})})],

{{\displaystyle u}}_j^{(\mathrm{3})}={{\displaystyle u}}_j^{(\mathrm{2})}+{\displaystyle \frac{\Delta t}{\mathrm{12}}}[-C({{\displaystyle u}}_j^{(\mathrm{0})})-C({{\displaystyle u}}_j^{(\mathrm{1})})+\mathrm{8}C({{\displaystyle u}}_j^{(\mathrm{2})})],

where

{{\displaystyle u}}_j^{(\mathrm{0})}:{{\displaystyle \bar{u}}}_j^n,{{\displaystyle u}}_j^{(\mathrm{3})}:={{\displaystyle \bar{u}}}_j^{n+\mathrm{1}}

, and

(57)C({{\displaystyle u}}_j^{(i)}):=-{\displaystyle \frac{{{\displaystyle F}}_{j+\frac{\mathrm{1}}{\mathrm{2}}}^{*}-{{\displaystyle F}}_{j-\frac{\mathrm{1}}{\mathrm{2}}}^{*}}{{{\displaystyle \Delta x}}_j}}+s(t^i),for\hspace{0.17em}i=\mathrm{0},\mathrm{1},\mathrm{2},

where s(t^{\mathrm{0}}):=s(t^n), s(t^{\mathrm{1}}):=s\left(t^n+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\Delta t\right) and s(t^{\mathrm{2}}):=s\left(t^n+{\displaystyle \frac{\mathrm{8}}{\mathrm{12}}}\Delta t\right).

For the numerical stability through time, it is necessary that the time step \Delta t satisfies the CFL (or the Courant) condition for the central schemes as:

(58){\displaystyle \frac{\Delta t}{{{\displaystyle \min }}_i({{\displaystyle \Delta x}}_i)}}{{\displaystyle \max }}_j\left|{{\displaystyle \lambda }}_j\right|\le {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},

where {{\displaystyle \lambda }}_j denotes the eigenvalues of the Jacobian matrix and {{\displaystyle \Delta x}}_i is length of the cell {{\displaystyle I}}_i. In practice, smaller values than {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}} is usually needed.

5.3 Complexity of computations

A direct solution (with the matrix inversion operator) of the linear system in Eq. (36) of size N\times N, where N=Q+n, requires O(N^{\mathrm{2}}) operations. Solution of this linear system is nearly the most frequent operation in RBF-based simulations: for M locations, it becomes O(M\hspace{0.17em}{N}^{\mathrm{2}}). Evaluations of s^j(x) in Eq. (32) at M points {{\displaystyle x}}_i need O(MN) operations; monotonization of s^j(x) by the nonlinear limiter (Eq. (39)) also increases this computation cost. High computation cost of RBFs were mentioned at first by Franke [92]; several approaches have then been developed to reduce the computational complexity [93‒95], for example in Ref. [96], inversion complexity is reduced from O(N^{\mathrm{2}}) to O(N\log (N)) operations and O(MN) of s^j(x) to O(M+N).

In our study, for each RBF, three successive cells (i.e., Q=\mathrm{3}) with zero order of polynomials (n=\mathrm{0}) are used for interpolations: that is N=\mathrm{3}. This means number of locations M control the computational complexity and in the worst case, the computational complexity is O({\mathrm{3}}^{\mathrm{3}}M) for the inversion operator and O(\mathrm{3}M) for s^j({{\displaystyle x}}_i). It should be mentioned that, in central high resolution schemes, values of reconstructed solutions at cell edges are needed (Eq. (55)); so, for {{\displaystyle M}}_c cells, M=\mathrm{2}{{\displaystyle M}}_c. This means that the computational complexity is directly controlled by cell numbers and the cost can be controlled by effective adaptation procedures. For this purpose, the fast MRA-based cell/grid adaptation is used in this work. For 1-D data of length M, wavelet transforms (with the pyramidal algorithm) use M operations and so they are fast and effective. In Section 7, by some numerical examples, it is shown that in the worst case, compression ratio of cells by the adaptation procedure is less than 15 percent, and this confirms effectiveness of adaptive solvers. In the following, we try to estimate computational complexity of the proposed central scheme by a numerical approach.

5.3.1 Complexity of the RBF-based central high resolution scheme

Following Ref. [97], we evaluate the complexity of the RBF-based central high resolution schemes by the numerical approach by some simulations. Let us assume the Burgers equation {{\displaystyle u}}_{,t}+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{u}^{\mathrm{2}}\right)=\mathrm{0},\hspace{0.17em}\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}x\in [\mathrm{0},\mathrm{1}], for x\in [\mathrm{0},\mathrm{1}] where u:=u(x,t); the boundary conditions (BCs) and initial condition (IC) are: BCs: u(\mathrm{0},t)=u(\mathrm{1},t)=\mathrm{0}, and IC: u(x,\mathrm{0})=\sin (\mathrm{2}\pi x)+{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\sin (\pi x). Simulations are done by the Gaussian-RBF with the optimum shape parameter, c=\mathrm{6.5}. The parameters used/provided in this study are:1) {{\displaystyle J}}_{\max } and {{\displaystyle J}}_{\min } denote the finest and coarsest resolution levels, respectively; 2) {{\displaystyle N}}_{{{\displaystyle J}}_{\max }} is the number of uniform cells in the finest resolution level and equal to {\mathrm{2}}^{{{\displaystyle J}}_{\max }}+\mathrm{1} (where the sampling step is \mathrm{1}/{\mathrm{2}}^{{{\displaystyle J}}_{\max }}); 3) {{\displaystyle N}}_{dt} is number of time-steps; 4) {{\displaystyle r}}_{N,\in }^{{{\displaystyle J}}_{\max }}:=\max ({{\displaystyle N}}_g)/{{\displaystyle N}}_{{{\displaystyle J}}_{\max }}denotes the minimum reduction (compression) factor in grid numbers during numerical simulations by the MRA-based grid adaptation method with the finest resolution level {{\displaystyle J}}_{\max } and the pre-defined threshold \in; and 5) {{\displaystyle C}}_{MS} denotes the computational time.

Complexities of the RBF-based central high resolution schemes without and with the adaptive solver are investigated separately in the following. It should be mentioned that since here three successive cells are used for the RBF formulation (i.e., Q=\mathrm{3}), closed-forms of both {{\displaystyle \overline{A}}}_{t,l}={{\displaystyle \mu }}_t^x{{\displaystyle \mu }}_l^y\phi (\Vert x-y\Vert ) and {{\displaystyle \mu }}_k^y\phi (\Vert x-y\Vert ) are derived at first and then used in all simulations.

i) Complexities of the RBF-based central high resolution schemes on uniform grids/cells

As mentioned above, as number of successive centers used in the RBFs is three (Q=\mathrm{2}q+\mathrm{1}=\mathrm{3}), the main source of computational complexity is number of locations M. In this study, different resolution levels j={{\displaystyle J}}_{\max } are used for simulations with location numbers

M={{\displaystyle N}}_{{{\displaystyle J}}_{\max }}={\mathrm{2}}^{{{\displaystyle J}}_{\max }}+\mathrm{1}

. To study only complexity of the central high resolution schemes, the time step dt is assumed to be constant for different uniform grids/cells of different cell-lengths. Results are measured and reported in Table 2 after 395 time steps at t=\mathrm{0.158}, when a right propagating shock wave developed for the first time. Theoretically, computational time {{\displaystyle C}}_{MS} decreases asymptotically by a factor {{\displaystyle r}}_c:={{\displaystyle ({{\displaystyle C}}_{MS})}}_{{{\displaystyle J}}_{\max }}/{{\displaystyle ({{\displaystyle C}}_{MS})}}_{{{\displaystyle J}}_{\max }-\mathrm{1}} . Regarding Table 2, this ratio is about 2 and so it is clear that the computational complexity asymptotically is O(M) operations.

ii) Complexities of the RBF-based central high resolution schemes on adapted grids/cells:

For numerical estimation of the computational complexity, the Burgers equation is re-considered. Assumptions in the simulations are: the Gaussian-RBF is used with C=\mathrm{6.5}; the wavelet threshold is \in =\mathrm{0.75}\times {\mathrm{10}}^{-\mathrm{2}}; the CFL value for all simulations is constant and equal to {\displaystyle \frac{\Delta t}{{{\displaystyle \min }}_i({{\displaystyle \Delta x}}_i)}}{{\displaystyle \max }}_i\left|{{\displaystyle \lambda }}_j\right|=\mathrm{0.25}; the coarsest level {{\displaystyle J}}_{\min } keeps constant and equal to 5 for all simulations (with the sampling step \Delta x=\mathrm{1}/{\mathrm{2}}^{{{\displaystyle J}}_{\min }}); to capture fine resolution effects, {{\displaystyle J}}_{\max } is increased one-by-one from 9 to 13; the MRA-based grid adaptation is done after each time step during time; closed-forms of both {{\displaystyle \mu }}_i^x{{\displaystyle \mu }}_l^y\phi and {{\displaystyle \mu }}_k^y are used in all simulations. Results are provided in Table 3 at t=\mathrm{0.158} (the time of developing the first discontinuity in solutions).

The relative number {{\displaystyle r}}_{N,\in }^{{{\displaystyle J}}_{\max }} decreases asymptotically by a factor {{\displaystyle r}}_j:={{\displaystyle r}}_{N,\in }^{{{\displaystyle J}}_{\max }-\mathrm{1}}/{{\displaystyle r}}_{N,\in }^{{{\displaystyle J}}_{\max }}; the results offer that it is about 2. The total number of adapted cells increases by a factor q:={\displaystyle \frac{{{\displaystyle [\max ({{\displaystyle N}}_g)]}}_{{{\displaystyle J}}_{\max }-\mathrm{1}}}{{{\displaystyle [\max ({{\displaystyle N}}_g)]}}_{{{\displaystyle J}}_{\max }}}}\approx {{\displaystyle \mathrm{2}r}}_j. The ratio {{\displaystyle r}}_j is around 0.5 and so {{\displaystyle q}}_l\approx \mathrm{1}. This means that the computational complexity in the spatial domain does not vary considerably by adding an extra finer resolution level for capturing fine responses. This reveals importance of adaptive solvers in the spatial-domain. On the other hand, by adding an additional refinement, due to the CFL condition, the time step decreases by a factor 0.5. Hence, the total computational complexity in the spatio-temporal domains would increase by {{\displaystyle \mathrm{2}q}}_l. This can be confirmed by the {{\displaystyle C}}_{MS} values as {{\displaystyle r}}_c:={{\displaystyle ({{\displaystyle C}}_{MS})}}_{{{\displaystyle J}}_{\max }}/{{\displaystyle ({{\displaystyle C}}_{MS})}}_{{{\displaystyle J}}_{\max }-\mathrm{1}} is around 2.

5.4 Computation costs of different third-order central high resolution schemes

Let reconsider the above-mentioned Burgers problem on uniform cells; the computational time of the RBF-based central high resolution scheme is presented in Table 2 for the Gaussian-RBF with c=\mathrm{6.5} after 395 time step at t=\mathrm{0.158}. Let us consider the third order central high resolution method of Liu and Tadmor (LT) [38] and the central-WENO (CWENO) scheme with the three-point stencil [40,41] as two polynomial-based third order central schemes. The LT scheme uses average interpolating parabolic polynomials monotonized by a nonlinear scaling limiter; the CWENO method employs a convex combination of two one-sided (left and right) and one central polynomials (with the average interpolating feature) to prevent developing of spurious oscillations. By these methods, computational times are presented in Table 4 after 395 time steps at t=\mathrm{0.158}. The results confirm that the RBF-based method is slower than both the LT and CWENO schemes on average by a factor of 2.55. However, it should be mentioned that the main advantage of the central RBF-based schemes is their stable formulations for handling of non-uniform cells with minimal spurious oscillations (in their numerical solutions).

5.5 Extension to 2-D problems

The 1-D RBF-based monotonized reconstruction and corresponding central high resolution scheme can be extended to higher dimensions. Depending of grid/cell configurations, there are some approaches:

1) Structured 2-D grids/cells: for this case, the dimension-by-dimension discretization can be used. Then the 1-D monotonized RBF-based interpolation can independently be performed for each dimension without any modification. The MRA-based grid/cell adaptation is also performed by the 1-D algorithm,

2) Scattered or unstructured 2-D grids/cells: the RBF-based reconstruction on triangular cells was successfully studied in Ref. [58]. For the monotonization step, however, it is necessary to update formulation of the (nonlinear) rescaling limiter in Eq. (39) to account all of the surrounding cells, see e.g., Ref. [98]. Also, MRA on scattered data was successfully developed and integrated for data analysis and PDE simulations, e.g., Ref. [99].

6 Implementation of RBF-based central high resolution schemes on MRA-based adapted grids

Let {{\displaystyle \bar{u}}}_j^n:=\bar{u}({{\displaystyle x}}_j,t^n) denote average solutions of a hyperbolic system on cells {{\displaystyle I}}_j. Solution {{\displaystyle \bar{u}}}_j^{n+\mathrm{1}} can be evaluated by the following procedure:

1) Determine the MRA-based adaptive grid from {{\displaystyle \bar{u}}}_j^n by using the D-D adaptive wavelet transform. Points without {{\displaystyle \bar{u}}}_j^n values are estimated by a local interpolation (e.g., by using splines),

2) Compute the RBF-based reconstruction over each cell and evaluate the piecewise interpolation s(x) (after choosing a good/optimum value of the shape parameter),

3) Monotonize s(x) by the nonlinear scaling limiter {{\displaystyle \theta }}_j and update s(x) as Eq. (43),

4) Evaluate {{\displaystyle u}}_{j\pm \mathrm{1}/\mathrm{2}}^L and {{\displaystyle u}}_{j\pm \mathrm{1}/\mathrm{2}}^R based on the updated s(x); then determine F({{\displaystyle u}}_{j\pm \mathrm{1}/\mathrm{2}}^L), and F({{\displaystyle u}}_{j\pm \mathrm{1}/\mathrm{2}}^R) and

{{\displaystyle a}}_{j\pm \mathrm{1}/\mathrm{2}}

,

5) Compute the reconstructed fluxes (in the semi-discrete formulation), i.e., {{\displaystyle F}}_{j\pm \mathrm{1}/\mathrm{2}}^{*},

6) Evaluate {{\displaystyle \bar{u}}}_j^{n+\mathrm{1}} by the third-order TVD Runge-Kutta scheme. In this stage, it is necessary to repeat steps 2 to 5 for two other internal steps; evaluate {{\displaystyle u}}_j^{(i)} and C({{\displaystyle u}}_j^{(i)}) at the end of each step, denoted by i. Using these values for i\{\mathrm{0},\mathrm{1},\mathrm{2}\}, calculate finally {{\displaystyle \bar{u}}}_j^{n+\mathrm{1}},

7) Set n:=n+\mathrm{1} and go back to step 1.

In practice, to have a cost-effective computation, the grid can be adapted after some time steps, depending on the velocity of moving fronts.

7 Numerical examples

The following examples are to study the effectiveness of the proposed method in the coupled thermoelasticity problems (1-D first order hyperbolic systems). Examples contain the GN and LS formulations with constant or variable coefficients.

The main assumptions are: 1- applying the D-D interpolating wavelet of order 3; 2- utilizing adapted grids with the finest and coarsest resolution levels: {{\displaystyle J}}_{\max }=\mathrm{11} and {{\displaystyle J}}_{\min }=\mathrm{5}; this means that the decomposing levels are six; 3- assuming that {{\displaystyle N}}_s={{\displaystyle N}}_c=\mathrm{2} in the MRA-based grid modification stage; 4- repeating the re-adaptation process at every time step; 5- integrating in time by the TVD Runge-Kutta solvers. For the RBF-based and the KT schemes, the Runge-Kutta schemes with second and third orders are used, respectively; 6- utilizing of the G-RBF in all simulations and MQ- and IMQ- RBFs in Example 1; 7- employing of three successive cells for RBF-based reconstructions (i.e.: q=\mathrm{1} in Eq. (31)); 8- using the generalized MINMOD flux limiter for KT solver where \widehat{\theta }=\mathrm{2} (see Appendix A).

7.1 Example 1: GN and LS one-dimensional generalized thermoelastic problems with constant coefficients

A domain is assumed with possible variations only in the x direction: one-dimensional problem. The assumptions are: 1) both the thermal source (r) and body force ({{\displaystyle f}}_i) are zero; 2) \mathrm{0}\le x\le \mathrm{2}; 3) the initial conditions (ICs) are zero, that is: v(x,t=\mathrm{0})=\mathrm{0}, \sigma (x,\mathrm{0})=\mathrm{0}, \theta (x,\mathrm{0})=\mathrm{0} and q(x,\mathrm{0})=\mathrm{0}; 4) the boundary conditions (BCs) are:{{\displaystyle v}}_{,x}(\mathrm{0},t)=\mathrm{0}, \theta (x=\mathrm{0},t)=H(t) (where H(t) denotes the Heaviside function), \sigma (\mathrm{0},t)=\mathrm{0} (a free BC) and {{\displaystyle q}}_{,x}(\mathrm{0},t)=\mathrm{0} (an isolated BC); at x=\mathrm{2}, the problem has fix and isolated BCs, as: v(\mathrm{2},t)=\mathrm{0}\theta (\mathrm{2},t)=\mathrm{0}, {{\displaystyle \sigma }}_{,x}(\mathrm{2},t)=\mathrm{0} and {{\displaystyle q}}_{,x}(\mathrm{0},t)=\mathrm{0}; 5) material parameters of the domain are: {{\displaystyle v}}_e=\mathrm{1}, {{\displaystyle v}}_t=\mathrm{2}, \rho =\mathrm{1}, \lambda =\mathrm{0.5}, \mu =\mathrm{0.25}, v=\mathrm{0.02}, a=\mathrm{0.05}, {{\displaystyle \tau }}_{\mathrm{0}}=\mathrm{0.5} and {{\displaystyle \alpha }}_t={\mathrm{10}}^{-\mathrm{2}}; 6) the reference temperature is {{\displaystyle T}}_{\mathrm{0}}=\mathrm{293.15} K (or 20°C); 7) the CFL condition is assumed to be: {\displaystyle \frac{{{\displaystyle \max }}_i\left|{{\displaystyle \bar{\lambda }}}_i\right|}{\Delta x/\Delta t}}=\mathrm{0.1}, where

\Delta x=\mathrm{1}/{\mathrm{2}}^{{{\displaystyle J}}_{\max }}

and {{\displaystyle \bar{\lambda }}}_i denote the eigenvalues of the Jacobian matrix A; 8) for grid adaptations, thresholds \in ={\mathrm{10}}^{-\mathrm{4}} and \in =\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{4}} are used for GN and LS problems, respectively.

For imposing the BCs, the concept of ghost cells is used. For the KT and RBF-based schemes respectively one and two ghost cells are used for each side. Let us denote the first and last cells on the original (non-extended) domain respectively by {{\displaystyle I}}_{\mathrm{1}} and {{\displaystyle I}}_N (and this means {{\displaystyle x}}_{\mathrm{1}-\mathrm{1}/\mathrm{2}}=\mathrm{0} and {{\displaystyle x}}_{N+\mathrm{1}/\mathrm{2}}=\mathrm{2}); for the state variables, based on BCs, the extended values with two ghost cells for each side are:

\mathrm{F}\mathrm{o}\mathrm{r}\hspace{0.17em}:{{\displaystyle I}}_{\mathrm{0}}:{{\displaystyle v}}_{\mathrm{0}}={{\displaystyle v}}_{\mathrm{1}},{{\displaystyle \sigma }}_{\mathrm{0}}=-{{\displaystyle \sigma }}_{\mathrm{1}},{{\displaystyle \theta }}_{\mathrm{0}}=\mathrm{2}H(t)-{{\displaystyle \theta }}_{\mathrm{1}},{{\displaystyle q}}_{\mathrm{0}}{{\displaystyle =q}}_{\mathrm{1}},

\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}:{{\displaystyle I}}_{-\mathrm{1}}:{{\displaystyle v}}_{-\mathrm{1}}={{\displaystyle v}}_{\mathrm{2}},{{\displaystyle \sigma }}_{-\mathrm{1}}=-{{\displaystyle \sigma }}_{\mathrm{2}},{{\displaystyle \theta }}_{-\mathrm{1}}=\mathrm{2}H(t)-{{\displaystyle \theta }}_{\mathrm{2}},{{\displaystyle q}}_{-\mathrm{1}}{{\displaystyle =q}}_{\mathrm{2}},

\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}:{{\displaystyle I}}_{N+\mathrm{1}}:{{\displaystyle v}}_{N+\mathrm{1}}={{\displaystyle v}}_N,{{\displaystyle \sigma }}_{N+\mathrm{1}}=-{{\displaystyle \sigma }}_N,{{\displaystyle \theta }}_{N+\mathrm{1}}=-{{\displaystyle \theta }}_N,{{\displaystyle q}}_{N+\mathrm{1}}{{\displaystyle =q}}_N,

(59)\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}:{{\displaystyle I}}_{N+\mathrm{2}}:{{\displaystyle v}}_{N+\mathrm{2}}={{\displaystyle v}}_{N-\mathrm{1}},{{\displaystyle \sigma }}_{N+\mathrm{2}}=-{{\displaystyle \sigma }}_{N-\mathrm{1}},{{\displaystyle \theta }}_{N+\mathrm{2}}=-{{\displaystyle \theta }}_{N-\mathrm{1}},{{\displaystyle q}}_{N+\mathrm{2}}{{\displaystyle =q}}_{N-\mathrm{1}}.

For the KT and RBF-based schemes corresponding numerical results for the medium with the GN description are presented in Fig. 2 at t=\mathrm{0.73} where the G-RBF is used with c=\mathrm{1}. In each illustration, distributions of adapted grids in different resolution levels resulted from different solvers are also presented. In Fig. 2, hollow shapes denote adaptive numerical results and the solid line is the exact [100] solution for the GN theory; also {\sigma }^{*}:=\sigma /\beta and v^{*}:={\displaystyle \frac{v\lambda +\mathrm{2}\mu }{\beta \hspace{0.17em}{{\displaystyle v}}_e}} are normalized state variables. In Fig. 2, variation of condition numbers (\kappa) is also presented; they increase rapidly on adapted cells/grids and this can increase numerical errors (this will be shown).

Performance of different RBFs around discontinuous solutions is presented in Fig. 3 for different c values at times \mathrm{0.031}, \mathrm{0.06} and \mathrm{0.12}. The results offer that: i) the G-RBF is not so sensitive to c values and adapted cells/grids; ii) small c values lead to localized spurious oscillations (for the G- and MQ-RBFs) or even instability (for the IMQ-RBF) due to large values of condition numbers, \kappa. So, for an optimum shape parameter of small value, using of larger values can be recommended. For this reason, for the IMQ-RBF with {{\displaystyle c}}_{opt}=\mathrm{4.5} (see Fig. 1), a larger value, i.e., c=\mathrm{6} is chosen as a good value. For different RBFs with their good values of c, numerical results and their condition numbers are presented in Fig. 4 at t=\mathrm{0.73}. The results provide that: 1) adaptive solvers are stable for all the RBFs; 2) adapted cells concentrate properly around high gradient zones; 3) around adapted cells, \kappa values increase considerably: the main source of numerical errors. These variation patterns reveal again importance of proper selection of c values.

Regarding the LS theory, adapted solutions with the RBF-based and KT solvers are provided in Fig. 5 where the G-RBF is used with c=\mathrm{1}; there, hollow shapes denote adaptive numerical results and the solid line is the asymptotic [24,25] solution for the LS theory where {\sigma }^{*}:=\sigma /\beta and v^{*}:={\displaystyle \frac{v\lambda +\mathrm{2}\mu }{\beta \hspace{0.17em}{{\displaystyle v}}_e}}. Variation of \kappa is also presented for the RBF-based solver; it is inversely proportional to cell sizes.

It should be mentioned that asymptotic solutions valid for small times; hence, difference between numerical and asymptotic solutions increase through time. For the LS problem, this is shown in Fig. 6 for the solution v^{*} at times 0.077, 0.154 and 0.232; in this figure, the solid lines and hollow shapes denote asymptotic and RBF-based numerical solutions, respectively.

Variations of number of adapted grid points, {{\displaystyle N}}_g, through time are illustrated in Fig. 7 for both the KT and the RBF-based central high resolution schemes (by the G-RBF with c=\mathrm{1}). In this figure, also, total variations (TV) of solution v, resulted from the RBF-based solvers, are presented; where the function TV of discrete data {{\displaystyle u}}_i is defined as:

(60)TV(u):={{\displaystyle \mathrm{\Sigma }}}_i\left|{{\displaystyle u}}_{i+\mathrm{1}}-{{\displaystyle u}}_i\right|,

Figure 7 offers that: i) variations of TV are bounded; ii) number of grid points in the finest resolution is 2049 points while number of adapted grids, {{\displaystyle N}}_gs are less than 300 for all simulations. This means that the compression rate is less than \mathrm{0.15} and this confirms effectiveness of the adaptive solvers.

Based on the above-mentioned comparisons between the KT and RBF-based central schemes, the adaptive RBF-based simulations lead to: i) smaller numerical dissipations and so more accurate results in comparison with the KT scheme; ii) solutions with possible localized spurious oscillations around discontinuities; iii) more mobilized adapted grid points especially of finer resolutions with respect to the KT method; iv) proper concentration of adapted grid points around high gradient zones and discontinuities; v) the TVB feature; vi) cost effective solvers.

7.2 Example 2:

In the previous example, effects of the high-frequency impulse BC are studied. In this example a sin-like BC is considered to simulate smooth but high-gradient solutions. So, it is assumed that:

(61)\theta (x=\mathrm{0},t)=H({{\displaystyle a}}_{\mathrm{0}}-t)sin\left({\displaystyle \frac{\pi t}{{{\displaystyle a}}_{\mathrm{0}}}}\right),

where {{\displaystyle a}}_{\mathrm{0}}=\mathrm{0.12}.

The thermoelastic material has the GN feature. Remaining BCs, ICs and parameters are the same as the previous example. Numerical results and corresponding adapted grid points are presented in Figs. 8 and 9 respectively for the KT and RBF-based schemes; the RBF is the Gaussian one with c=\mathrm{1}. For the RBF-solvers, also, variation of \kappa (condition number) values is provided. In these figures, the solid lines and hollow shapes are exact [100] and numerical solutions, respectively. The results offer that both schemes can properly handle high gradient smooth solutions; however, the RBF-based scheme is more accurate and can mobilize finer resolutions. These fine resolutions, however, can increase considerably the condition numbers and so using of a proper value of c is important.

7.3 Example 3:

Let us reconsider example 1, but in this case with variable coefficients. That is \kappa, {\kappa }^{*} and {{\displaystyle c}}_pdepend linearly to the temperature increment parameter, \theta. Considering Eq. (13), it is assumed: {{\displaystyle \chi }}_{\mathrm{1}}={{\displaystyle \chi }}_{\mathrm{2}}=\mathrm{0.25}. For numerical simulations, the G-RBF is used with c=\mathrm{1} and the remaining parameters are the same as example 1.

Numerical results and corresponding adapted grids for the GN and LS problems are illustrated in Figs. 10 and 11, respectively. There, the hollow shapes and solid lines are numerical and asymptotic [24,25] solutions, respectively. In both figures, variations of \kappa (condition number) values of RBF-solvers are also provided.

Variations of {{\displaystyle N}}_g through time for both the KT and RBF-based solvers (the Gaussian one with c=\mathrm{1}) along with TVs of v for the RBF-based scheme are presented in Fig. 12.

Numerical results again confirm that the RBF-scheme leads to more accurate results (with some localized non-physical oscillations), more mobilized adapted points and the TVB feature. Adapted points also locate properly around high gradient regions. Please notice that since asymptotic solutions are used for comparisons, these solutions have some differences with the numerical results (see v^{*} in Fig. 11).

8 Conclusions

In this study, a higher order central high resolution scheme is introduced by using of RBFs in the reconstruction stage. The method has the TVB feature and this is because: i) over each cells, RBFs have optimum recovery with minimized roughness; ii) at cell edges, the monotonicity is checked/imposed by using the scaling nonlinear limiter. Due to stability of RBFs on non-uniform cells/grids, the proposed central scheme is also stable on such grids. The relationship between reconstructions with the RBFs and polynomials is investigated by the Taylor series. And it is shown that the RBF-based method is a higher order scheme (at least) of third order accuracy.

This RBF-based central scheme is used for simulation of some generalized thermoelasticity problems with constant and variable coefficients. The GN and LS theories are presented as a unified formulation. For systems with constant and variable coefficients, corresponding coupled first order hyperbolic systems, Jacobian matrices and eigenvalues are derived. This is because, the maximum absolute value of eigenvalues (regardless of its propagation direction) is essential in formulations of the higher order central high resolution schemes.

Computational complexities are provided and also computational cost of the central RBF-based solver is compared with two other polynomial-based third-order high resolution schemes. These studies reveal that the main drawback of higher-order central high resolution schemes especially the RBF-based solvers is that they are expensive. In this regard, developing of adaptive high-order solvers is important. In this study, the MRA-based approach is used for cell/grid adaptations. The adaptations are done by the D-D interpolating wavelets. On the other hand, by using finer resolutions, condition numbers of RBF-based solvers increase considerably and this leads to numerical errors or even instabilities. Shape parameters of the RBFs can control trade-off between stability and accuracy. Hence, proper selection of the shape parameters is important. In this study, good/optimum values of the shape parameters for different RBFs are simply selected by some test functions and an error definition containing reconstruction errors at cell-edges. This numerical approach is done for the G, MQ and IMQ RBFs. It is shown that the MQ and IMQ RBFs are very sensitive for small values of the shape parameter (c) and cells of small lengths (\Delta x), and so a proper selection of c values is important.

To compare performance of TVB with TVD solvers (both on adapted grids), the second order KT central scheme is also used in numerical simulations of some benchmarks. For the adaptive KT central high resolution scheme see Appendix A. Both the TVB and TVD approaches are successfully integrated with the MRA-based grid adaptation. The numerical results confirm that:

1) Adapted points concentrate properly over high gradient zones and around discontinuities,

2) The RBF-based central scheme has the TVB property.

3) The RBF-based scheme mobilizes more adapted grids due to its higher-order accuracy and TVB features. The TVB property let localized spurious oscillations develop in numerical solutions,

4) The RBF-based scheme has more accuracy and less numerical dissipation than the KT method. It should be mentioned that the TVD criterion leads to the first order accuracy around discontinuities, while the TVB one yields high order accuracy everywhere in solutions ever in the vicinity of discontinuities,

5) Adaptive solutions are effective since less than 15 percent of grid points in the finest resolution level is mobilized in numerical simulations,

6) Thermoelastic coupled shock waves with finite velocities are properly captured for both the GN and LS formulations with constant or variable coefficients. Simulations are done without adding extra artificial viscosity (damping) around discontinuities or high-frequency impulse load [101],

7) Both impulsive and smooth-high-gradient BCs are successfully simulated by the RBF-based central scheme. In conclusion, this method is more accurate than the KT scheme.

References

[1] Chandrasekharaiah D S. Thermoelasticity with second sound: a review. Applied Mechanics Reviews, 1986, 39(3): 355–376

[2] Mallik S H, Kanoria M. Generalized thermoelastic functionally graded solid with a periodically varying heat source. International Journal of Solids and Structures, 2007, 44(22‒23): 7633–7645

[3] Tamma K K, Namburu R R. Computational approaches with applications to non-classical and classical thermomechanical problems. Applied Mechanics Reviews, 1997, 50(9): 514–551

[4] Mitra K, Kumar S, Vedevarz A, Moallemi M. Experimental evidence of hyperbolic heat conduction in processed meat. Journal of Heat Transfer, 1995, 117(3): 568–573

[5] Lord H W, Shulman Y. A generalized dynamical theory of thermoelasticity. Journal of the Mechanics and Physics of Solids, 1967, 15(5): 299–309

[6] Green A, Lindsay K. Thermoelasticity. Journal of Elasticity, 1972, 2(1): 1–7

[7] Green A, Naghdi P. Thermoelasticity without energy dissipation. Journal of Elasticity, 1993, 31(3): 189–208

[8] Green A E, Naghdi P M. A re-examination of the basic postulates of thermomechanics. Proceedings of the Royal Society of London. Series A, 1991, 432(1885): 171–194

[9] Green A, Naghdi P. On undamped heat waves in an elastic solid. Journal of Thermal Stresses, 1992, 15(2): 253–264

[10] Das B. Problems and Solutions in Thermoelasticity and Magneto-thermoelasticity. Springer, 2017

[11] Povstenko Y. Fractional Cattaneo-type equations and generalized thermoelasticity. Journal of Thermal Stresses, 2011, 34(2): 97–114

[12] Povstenko Y. Fractional thermoelasticity. In: Hetnarski R B, ed. Encyclopedia of Thermal Stresses. Springer, 2014, 1778–1787

[13] Ezzat M A, El-Karamany A S, Samaan A A. The dependence of the modulus of elasticity on reference temperature in generalized thermoelasticity with thermal relaxation. Applied Mathematics and Computation, 2004, 147(1): 169–189

[14] Youssef H M. Dependence of modulus of elasticity and thermal conductivity on reference temperature in generalized thermoelasticity for an infinite material with a spherical cavity. Applied Mathematics and Mechanics, 2005, 26(4): 470–475

[15] Aouadi M. Generalized thermo-piezoelectric problems with temperature-dependent properties. International Journal of Solids and Structures, 2006, 43(21): 6347–6358

[16] Othman M I, Kumar R. Reflection of magneto-thermoelasticity waves with temperature dependent properties in generalized thermoelasticity. International Communications in Heat and Mass Transfer, 2009, 36(5): 513–520

[17] Allam M N, Elsibai K A, Abouelregal A E. Magneto-thermoelasticity for an infinite body with a spherical cavity and variable material properties without energy dissipation. International Journal of Solids and Structures, 2010, 47(20): 2631–2638

[18] Abbas I A. Eigenvalue approach in a three-dimensional generalized thermoelastic interactions with temperature-dependent material properties. Computers & Mathematics with Applications (Oxford, England), 2014, 68(12): 2036–2056

[19] Xiong Q L, Tian X G. Transient magneto-thermoelastic response for a semi-infinite body with voids and variable material properties during thermal shock. International Journal of Applied Mechanics, 2011, 3(4): 881–902

[20] He T, Shi S. Effect of temperature-dependent properties on thermoelastic problems with thermal relaxations. Chinese Journal of Solid Mechanics, 2014, 27: 412–419

[21] Sherief H, Abd El-Latief A M. Effect of variable thermal conductivity on a half-space under the fractional order theory of thermoelasticity. International Journal of Mechanical Sciences, 2013, 74: 185–189

[22] Wang Y, Xue J. Asymptotic analysis of thermoelastic response in a functionally graded solid based on LS theory. International Journal of Material Science, 2016, 6(1): 35–40

[23] Wang Y, Liu D, Wang Q, Zhou J. Thermoelastic behavior of elastic media with temperature-dependent properties under transient thermal shock. Journal of Thermal Stresses, 2016, 39(4): 460–473

[24] Wang Y, Liu D, Wang Q, Zhou J. Asymptotic solutions for generalized thermoelasticity with variable thermal material properties. Archives of Mechanics, 2016, 68: 181–202

[25] Liang W, Huang S, Tan W, Wang Y. Asymptotic approach to transient thermal shock problem with variable material properties. Mechanics of Advanced Materials and Structures, 2017: 1–9

[26] Youssef H, El-Bary A. Thermal shock problem of a generalized thermoelastic layered composite material with variable thermal conductivity. Math Probl Eng, 2006, 2006: 1–14

[27] Yousefi H, Noorzad A, Farjoodi J, Vahidi M. Multiresolution-based adaptive simulation of wave equation. Applied Mathematics & Information Sciences, 2012, 6: 47S–58S

[28] Yousefi H, Ghorashi S S, Rabczuk T. Directly simulation of second order hyperbolic systems in second order form via the regularization concept. Communications in Computational Physics, 2016, 20(01): 86–135

[29] Latifi M, Kharazi M, Ovesy H. Nonlinear dynamic response of symmetric laminated composite beams under combined in-plane and lateral loadings using full layerwise theory. Thin-walled Structures, 2016, 104: 62–70

[30] Latifi M, Farhatnia F, Kadkhodaei M. Buckling analysis of rectangular functionally graded plates under various edge conditions using Fourier series expansion. European Journal of Mechanics. A, Solids, 2013, 41: 16–27

[31] Latifi M, Kharazi M, Ovesy H. Effect of integral viscoelastic core on the nonlinear dynamic behaviour of composite sandwich beams with rectangular cross sections. International Journal of Mechanical Sciences, 2017, 123: 141–150

[32] Yousefi H, Noorzad A, Farjoodi J. Simulating 2D waves propagation in elastic solid media using wavelet based adaptive method. Journal of Scientific Computing, 2010, 42(3): 404–425

[33] Godunov S K. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. Matematicheskii Sbornik, 1959, 89: 271–306

[34] Harten A, Engquist B, Osher S, Chakravarthy S R. Uniformly high order accurate essentially non-oscillatory schemes, III. Journal of Computational Physics, 1987, 71(2): 231–303

[35] Liu X D, Osher S, Chan T. Weighted essentially non-oscillatory schemes. Journal of Computational Physics, 1994, 115(1): 200–212

[36] LeVeque R J. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, 2002

[37] Kurganov A, Tadmor E. New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations. Journal of Computational Physics, 2000, 160(1): 241–282

[38] Liu X D, Tadmor E. Third order nonoscillatory central scheme for hyperbolic conservation laws. Numer Math, 1998, 79(3): 397–425

[39] Van Leer B. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method. Journal of Computational Physics, 1979, 32(1): 101–136

[40] Levy D, Puppo G, Russo G. Central WENO schemes for hyperbolic systems of conservation laws. Modélisation Mathématique et Analyse Numérique, 1999, 33(3): 547–571

[41] Levy D, Puppo G, Russo G. Compact central WENO schemes for multidimensional conservation laws. SIAM Journal on Scientific Computing, 2000, 22(2): 656–672

[42] Dehghan M, Jazlanian R. On the total variation of a third-order semi-discrete central scheme for 1D conservation laws. Journal of Vibration and Control, 2011, 17(9): 1348–1358

[43] Levy A, Ben-Dor G, Sorek S. Numerical investigation of the propagation of shock waves in rigid porous materials: development of the computer code and comparison with experimental results. Journal of Fluid Mechanics, 1996, 324(-1): 163–179

[44] Levy A, Ben-Dor G, Sorek S. Numerical investigation of the propagation of shock waves in rigid porous materials: flow field behavior and parametric study. Shock Waves, 1998, 8(3): 127–137

[45] Heuzé T. Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids. Journal of Computational Physics, 2017, 346: 369–388

[46] Berezovski A, Maugin G. Simulation of thermoelastic wave propagation by means of a composite wave-propagation algorithm. Journal of Computational Physics, 2001, 168(1): 249–264

[47] Berezovski A, Maugin G. Application of wave-propagation algorithm to two-dimensional thermoelastic wave propagation in inhomogeneous media. In: Toro E F, ed. Godunov Methods: Theory and Applications. Boston: Springer Science & Business Media, 2001, 109–116

[48] Dehghan M, Abbaszadeh M. The space-splitting idea combined with local radial basis function meshless approach to simulate conservation laws equations. Alexandria Eng J, 2017,

[49] Shokri A, Dehghan M. A meshless method using radial basis functions for the numerical solution of two-dimensional complex Ginzburg-Landau equation. Computer Modeling in Engineering & Sciences, 2012, 84: 333‒358

[50] Guo J, Jung J H. Radial basis function ENO and WENO finite difference methods based on the optimization of shape parameters. Journal of Scientific Computing, 2017, 70(2): 551–575

[51] Dehghan M, Shokri A. A numerical method for solution of the two-dimensional sine-Gordon equation using the radial basis functions. Mathematics and Computers in Simulation, 2008, 79(3): 700–715

[52] Bigoni C, Hesthaven J S. Adaptive WENO methods based on radial basis function reconstruction. Journal of Scientific Computing, 2017, 72(3): 986–1020

[53] Sonar T. Optimal recovery using thin plate splines in finite volume methods for the numerical solution of hyperbolic conservation laws. IMA Journal of Numerical Analysis, 1996, 16(4): 549–581

[54] Golomb M, Weinberger H F. Optimal approximations and error bounds. In: Langer R E, ed. On Numerical Approximation. Madison: The University of Wisconsin Press, 1958, 117–190

[55] Micchelli C A, Rivlin T J. A survey of optimal recovery. In: Micchelli C A, Rivlin T J, eds. Optimal Estimation in Approximation Theory. Springer, 1977, 1–54

[56] Hickernell F J, Hon Y. Radial basis function approximations as smoothing splines. Applied Mathematics and Computation, 1999, 102(1): 1–24

[57] Iske A, Sonar T. On the structure of function spaces in optimal recovery of point functionals for ENO-schemes by radial basis functions. Numer Math, 1996, 74(2): 177–201

[58] Aboiyar T, Georgoulis E H, Iske A. High order WENO finite volume schemes using polyharmonic spline reconstruction. In: Agratini O, Blaga P, eds. International Conference on Numerical Analysis and Approximation Theory. Cluj-Napoca: Babeş – Bolyai University, CASA CĂRŢII DE ŞTIINŢĂ, 2006, 113–126

[59] Guo J, Jung J H. A RBF-WENO finite volume method for hyperbolic conservation laws with the monotone polynomial interpolation method. Applied Numerical Mathematics, 2017, 112: 27–50

[60] Schaback R. Error estimates and condition numbers for radial basis function interpolation. Advances in Computational Mathematics, 1995, 3(3): 251–264

[61] Mallat S. A Wavelet Tour of Signal Processing. Academic press, 1999

[62] Alves M, Cruz P, Mendes A, Magalhaes F, Pinho F, Oliveira P. Adaptive multiresolution approach for solution of hyperbolic PDEs. Computer Methods in Applied Mechanics and Engineering, 2002, 191(36): 3909–3928

[63] Santos J, Cruz P, Alves M, Oliveira P, Magalhães F, Mendes A. Adaptive multiresolution approach for two-dimensional PDEs. Computer Methods in Applied Mechanics and Engineering, 2004, 193(3‒5): 405–425

[64] Cohen A, Kaber S, Müller S, Postel M. Fully adaptive multiresolution finite volume schemes for conservation laws. Mathematics of Computation, 2003, 72(241): 183–225

[65] Dahmen W, Gottschlich–Müller B, Müller S. Multiresolution schemes for conservation laws. Numer Math, 2001, 88(3): 399–443

[66] Harten A. Multiresolution algorithms for the numerical solution of hyperbolic conservation laws. Communications on Pure and Applied Mathematics, 1995, 48(12): 1305–1342

[67] Gottschlich-Miiller B, Miiller S. Application of multiscale techniques to hyperbolic conservation laws. In: Chen Z, Li Y, Micchelli C, Xu Y, eds. Advances in Computational Mathematics, Lecture Notes in Pure & Applied Mathematics. Gaungzhou: Marcel Dekker, Inc., 1998, 113–138

[68] Berres S, Burger R, Kozakevicius A. Numerical approximation of oscillatory solutions of hyperbolic-elliptic systems of conservation laws by multiresolution schemes. Advances in Applied Mathematics and Mechanics, 2009, 1: 581–614

[69] Holmström M. Solving hyperbolic PDEs using interpolating wavelets. SIAM Journal on Scientific Computing, 1999, 21(2): 405–420

[70] Yousefi H, Noorzad A, Farjoodi J. Multiresolution based adaptive schemes for second order hyperbolic PDEs in elastodynamic problems. Applied Mathematical Modelling, 2013, 37(12‒13): 7095–7127

[71] Donoho D L, Johnstone J M. Ideal spatial adaptation by wavelet shrinkage. Biometrika, 1994, 81(3): 425–455

[72] Wang J, Liu G. Radial point interpolation method for elastoplastic problems. In. ICSSD 2000: 1 st Structural Conference on Structural Stability and Dynamics. 2000, 703–708

[73] Wang J, Liu G. On the optimal shape parameters of radial basis functions used for 2-D meshless methods. Computer Methods in Applied Mechanics and Engineering, 2002, 191(23‒24): 2611–2630

[74] Fasshauer G E. Meshfree Approximation Methods with Matlab. World Scientific Publishing Co Inc., 2007

[75] Liu G R. Gu Y-T. An Introduction to Meshfree Methods and Their Programming. Springer Science & Business Media, 2005

[76] Wendland H. Scattered Data Approximation. Cambridge University Press, 2004

[77] Driscoll T A, Fornberg B. Interpolation in the limit of increasingly flat radial basis functions. Computers & Mathematics with Applications (Oxford, England), 2002, 43(3‒5): 413–422

[78] Fornberg B, Larsson E, Flyer N. Stable computations with Gaussian radial basis functions. SIAM Journal on Scientific Computing, 2011, 33(2): 869–892

[79] Wright G B, Fornberg B. Stable computations with flat radial basis functions using vector-valued rational approximations. Journal of Computational Physics, 2017, 331: 137–156

[80] Fornberg B, Wright G. Stable computation of multiquadric interpolants for all values of the shape parameter. Computers & Mathematics with Applications (Oxford, England), 2004, 48(5‒6): 853–867

[81] Fasshauer G E, Zhang J G. Scattered data approximation of noisy data via iterated moving least squares. Curves and Surfaces: Avignon, 2006

[82] Fasshauer G E, Zhang J G. On choosing “optimal” shape parameters for RBF approximation. Numer Algor, 2007, 45(1‒4): 345–368

[83] Aboiyar T, Georgoulis E H, Iske A. Adaptive ADER methods using kernel-based polyharmonic spline WENO reconstruction. SIAM Journal on Scientific Computing, 2010, 32(6): 3251–3277

[84] Kansa E, Carlson R. Improved accuracy of multiquadric interpolation using variable shape parameters. Computers & Mathematics with Applications (Oxford, England), 1992, 24(12): 99–120

[85] Fornberg B, Zuev J. The Runge phenomenon and spatially variable shape parameters in RBF interpolation. Computers & Mathematics with Applications (Oxford, England), 2007, 54(3): 379–398

[86] Kansa E J. Multiquadrics—A scattered data approximation scheme with applications to computational fluid-dynamics—II solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers & Mathematics with Applications (Oxford, England), 1990, 19(8‒9): 147–161

[87] Hardy R L. Multiquadric equations of topography and other irregular surfaces. Journal of Geophysical Research, 1971, 76(8): 1905–1915

[88] Rippa S. An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Advances in Computational Mathematics, 1999, 11(2/3): 193–210

[89] Sanyasiraju Y, Satyanarayana C. On optimization of the RBF shape parameter in a grid-free local scheme for convection dominated problems over non-uniform centers. Applied Mathematical Modelling, 2013, 37(12‒13): 7245–7272

[90] Fjordholm U S, Ray D. A sign preserving WENO reconstruction method. Journal of Scientific Computing, 2016, 68(1): 42–63

[91] Gottlieb S, Shu C W, Tadmor E. Strong stability-preserving high-order time discretization methods. SIAM Review, 2001, 43(1): 89–112

[92] Franke R. Scattered data interpolation: tests of some methods. Mathematics of Computation, 1982, 38: 181–200

[93] Powell M J D. Tabulation of thin plate splines on a very fine two-dimensional grid. In: Numerical Methods in Approximation Theory, Vol. 9. Springer, 1992, 221–244

[94] Arad N, Dyn N, Reisfeld D, Yeshurun Y. Image warping by radial basis functions: application to facial expressions. Graphical Models, 1994, 56(2): 161–172

[95] Powell M. Truncated Laurent expansions for the fast evaluation of thin plate splines. Numer Algor, 1993, 5(2): 99–120

[96] Roussos G, Baxter B J. Rapid evaluation of radial basis functions. Journal of Computational and Applied Mathematics, 2005, 180(1): 51–70

[97] Müller S. Adaptive Multiscale Schemes for Conservation Laws. Springer Science & Business Media, 2012

[98] Cueto-Felgueroso L, Colominas I. High-order finite volume methods and multiresolution reproducing kernels. Archives of Computational Methods in Engineering, 2008, 15(2): 185–228

[99] Iske A. Multiresolution Methods in Scattered Data Modelling. Springer Science & Business Media, 2004

[100] Abd El-Latief A M, Khader S E. Exact solution of thermoelastic problem for a one-dimensional bar without energy dissipation. ISRN Mech Eng, 2014, 2014: 1–6

[101] Guo P, Wu W H, Wu Z G. A time discontinuous Galerkin finite element method for generalized thermo-elastic wave analysis, considering non-Fourier effects. Acta Mechanica, 2014, 225(1): 299–307

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature