Inverse differential quadrature method: mathematical formulation and error analysis

Engineering systems are typically governed by systems of high-order differential equations which require efficient numerical methods to provide reliable solutions, subject to imposed constraints. The conventional approach by direct approximation of system variables can potentially incur considerable error due to high sensitivity of high-order numerical differentiation to noise, thus necessitating improved techniques which can better satisfy the requirements of numerical accuracy desirable in solution of high-order systems. To this end, a novel inverse differential quadrature method (iDQM) is proposed for approximation of engineering systems. A detailed formulation of iDQM based on integration and DQM inversion is developed separately for approximation of arbitrary low-order functions from higher derivatives. Error formulation is further developed to evaluate the performance of the proposed method, whereas the accuracy through convergence, robustness and numerical stability is presented through articulation of two unique concepts of the iDQM scheme, known as Mixed iDQM and Full iDQM. By benchmarking iDQM solutions of high-order differential equations of linear and nonlinear systems drawn from heat transfer and mechanics problems against exact and DQM solutions, it is demonstrated that iDQM approximation is robust to furnish accurate solutions without losing computational efficiency, and offer superior numerical stability over DQM solutions.

Engineering systems are typically governed by systems of high-order differential equations which require efficient numerical methods to provide reliable solutions, subject to imposed constraints. The conventional approach by direct approximation of system variables can potentially incur considerable error due to high sensitivity of high-order numerical differentiation to noise, thus necessitating improved techniques which can better satisfy the requirements of numerical accuracy desirable in solution of high-order systems. To this end, a novel inverse differential quadrature method (iDQM) is proposed for approximation of engineering systems. A detailed formulation of iDQM based on integration and DQM inversion is developed separately for approximation of arbitrary low-order functions from higher derivatives. Error formulation is further developed to evaluate the performance of the proposed method, whereas the accuracy through convergence, robustness and numerical stability is presented through articulation of two unique concepts of the iDQM scheme, known as Mixed iDQM and Full iDQM. By benchmarking iDQM solutions of high-order differential equations of linear and nonlinear systems drawn from heat transfer and mechanics problems against exact and DQM solutions, it is demonstrated that iDQM approximation is robust to furnish accurate solutions without losing computational efficiency, and offer superior numerical stability over DQM solutions.

Introduction
Engineering systems are typically governed by complex high-order differential equations which require numerical methods to provide accurate solutions. To approximate such systems, the domain of interest is discretized by means of interpolation (shape) functions and its higher derivatives, which are defined over a subdomain of interest called elements. Examples of methods available in this context are element-based methods such as finite-element method [1][2][3], boundary element method [4,5], finite difference method [6][7][8] or finite volume method [9,10]. On the other hand, a class of high-order mesh-free methods such as radial basis function networks (RBFN) [11][12][13], element-free Galerkin method [14][15][16], diffuse element method [17], or highorder collocation methods such as the Chebyshev method [18][19][20] and the differential quadrature method (DQM) [21][22][23][24][25][26][27], have been widely applied for solving engineering and science problems. In all of these methods, target derivatives of an arbitrary function are directly approximated as a weighted sum of the function in the domain. According to Mai-Duy & Tran-Cong [28], in the process of differentiation, errors of function approximation may amplify significantly as influenced by local effects of the approximant. With respect to analysis of engineering structures, such error amplification can affect the accuracy of high-order secondary variables including strains, stresses, moments, and shear forces. To this effect, indirect radial basis function networks (IRBFN) were proposed in [29]. In contrast to the RBFN approach, IRBFN approximate a derivative function using RBFN and then recover the original function by integration, which is less sensitive to noise.
DQM as proposed by Bellman and others [21,22] has received widespread attention in the research community due to spectral accuracy and fast convergence that make it desirable for vast engineering applications. For example, in the field of structural mechanics, the DQM approach has been explored for static analysis [30][31][32], free vibration analysis [33][34][35], buckling analysis [29,36] and large deflection post-buckling analysis [37]. In the context of fluid mechanics applications, DQM has been widely applied to obtain solutions of convection problems [38][39][40] and Navier-Stokes equations [41,42], heat transfer analysis [43,44], and magnetohydrodynamic duct flow problems [45]. It is instructive to note that in [38][39][40]42] the so-called localized DQM is adopted in combination with RBF since this strategy allows versatility in using meshes or meshless systems [46]. Financial engineering is another area where the merits of DQM have been explored for computational analysis [47].
To generalize DQM for engineering applications, Shu [23] proposed a general approach which admits the use of different base polynomials for computation of DQM weights. According to Shu's approach [23], computation of the differential quadrature (DQ) weights for firstorder and higher-order derivatives can be generalized by appropriate choice of base vectors in the linear vector space. Furthermore, Quan & Chang [48] as well as Wang [49] contended that explicit formulation of the weighting coefficients by using test functions can overcome the ill-conditioning arising from a large number of grid points. In this regard, to allow for explicit computation of the DQ weights, setting the base vectors as base polynomials (such as Lagrange polynomials, Chebyshev polynomials or Legendre polynomials) in a linear vector space in the so-called polynomial-based differential quadrature (PDQ) method is recommended [23,48,49]. To mitigate the sensitivity of DQM solutions to grid spacing to the Runge phenomenon, Wang [49] noted that non-uniform grid distribution is required to obtain reliable solutions. Readers are referred to Wang [49] for classical examples of non-uniform grid distributions that guarantee good accuracy, numerical stability and fast convergence. On the other hand, for applications like wave propagation in space where uniformly distributed grid points are required for accurate numerical estimation, the localized DQM (LDQM), which applies DQM approximations within a small neighbourhood of the point of interest, is crucial to keep the balance between the accuracy and stability of the numerical estimates [25]. In recent times, remarkably in computational fluid mechanics applications, the so-called RBF-based DQ method, as well as the LDQM variant, are increasingly being adopted to extend the merits of the DQM, as already noted in [38][39][40]42]. where c k represents constant weights to be determined and ξ k are the linearly independent basis vectors in the N-dimensional linear vector space, V N . Adopting Shu's general approach [23], some sets of base polynomials can be selected to determine the weights, c k . Furthermore, according to Shu [23], the numerical difficulty of determining the weights as N becomes large can be eliminated by considering Lagrange interpolation polynomials as basis vectors in which P n (ξ ) is then approximated by where E(ξ ) is the polynomial approximation error, and l i (ξ ) is a (N − 1) degree Lagrange polynomial defined as , and f (ξ i ) is the functional value at a discrete point i. In line with DQM routines outlined in [23], the first derivative of the function f (ξ ), f (1) i , at a discrete point i in the domain ξ ∈ [a b] can be realized as follows, ij , for i = j. In equation (1.4), E dif(1) (ξ ) is the approximation error due to first-order numerical differentiation of f (ξ ), which is given in [23] as where K is a constant implied from Shu [23]. For an arbitrary derivative of order m, the DQM approximation is characterized by a weighting coefficient a ij which is evaluated based on the recursive formula, In a compact form, the DQM approximation of mth-order derivative of f (ξ ) is represented thus to equation (1.4), and E dif(m) is the approximation error of the mth-order numerical differentiation of F given as Typically, E dif(m) << F (m) , so equation (1.8) can be rewritten as where c (m) 0 is the mean approximation error of the mth-order numerical differentiation expressed as (1.10)

Inverse differential quadrature method
This section introduces the mathematical formulation for the iDQM. The idea of iDQM involves functional approximation of high-order derivatives of f (ξ ) and the subsequent recovering of loworder derivatives by integration of the high-order derivatives. Different approaches to achieve this aim in a computationally efficient and numerically stable manner are presented in the following.
(a) First-order inverse differential quadrature method-by-integration Suppose we let the first derivative of a function f (ξ ), f (1) (ξ ), be approximated for a fixed N as in equation (1.2), is the first derivative of f (ξ ) at a discrete point and E 1 (ξ ) is the polynomial approximation error for f (1) . The original function f (ξ ) can then be recovered from f (1) by integrating equation (2.1) to get where H i (ξ ) is a Nth-order polynomial function, c 0 is the constant of integration and E int(1) 1 (ξ ) = E 1 (ξ ) dξ . In a compact form, equation (2.2) gives It is convenient to recast equation (2.3) as: (b) Higher-order inverse differential quadrature method-by-integration Let the mth-order derivative of a continuous function f (ξ ) be f (m) (ξ ), which is approximated as in equation (1.2), is the polynomial approximation error for f (m) . To obtain f (m) (ξ ) from equation (2.5) requires mth order integration, which leads to Note: the subscript and superscript of E int(m) m refers, respectively, to the order of function, i.e. f (m) , and the order of integration operation.
In consistency with equations (2.2) and (2.3), equation (2.6) can be recast in a compact form as (c) Computational aspects of inverse differential quadrature method-by-integration The matrix of integral coefficients, H (m) , can be computed by analytical integration, which requires a costly symbolic computational operation in MATLAB. Besides, analytical integration may not be feasible for some vector basis functions, making this approach cumbersome and undesirable. On the other hand, Gaussian integration can be used to compute H (m) efficiently. However, the numerical accuracy and stability of this operation depends on the number of Gauss points, which is typically chosen intuitively. The additional computational variable manifested by Gaussian integration increases the complexity of the iDQM-by-integration especially for highorder functions. On this basis, to compute H (m) , we seek a novel alternative that, respectively, resolves the inefficiency and numerical instability bottlenecks of analytical and Gaussian integrations yet sufficiently preserves the accuracy of both methods.
(d) First-order inverse differential quadrature method-by-inversion To avoid the deficiency imposed by direct integration or Gaussian integration in the computation of H (m) , we now describe an alternative formulation which is computationally efficient and numerically stable to compute the iDQM coefficient matrix. Assuming the coefficient matrix, D, is invertible, equation (1.7) can be rearranged to make the vector function, F, the subject, in case of first order as (2.10) or alternatively as in equation (1.9), where D (1) = D −1 and c 0 is the is the mean error distribution for first-order iDQM-by-inversion, which is expressed as To compare equation ( Now, evaluating f (ξ p ) using equations (2.3) and (2.11), the following relation can be established, where [D (1) ] p ∈ R 1×N is the pth row vector of matrix D (1) . Rearranging equation (2.14) in terms of c 0 gives which is recast as where . Equation (2.17) is the equivalent of equation (2.4), which is numerically stable and computationally efficient for the approximation of f (ξ ) from its first derivative.
(e) Higher-order inverse differential quadrature method-by-inversion To obtain higher-order iDQM formulae from the DQM counterpart, we revisit equation (2.17) for the first-order derivative, explicitly expressed as Analogously, we can write derivations of first and second derivatives of f (ξ ) in terms of the second and third derivatives, respectively, as: and are error estimates implied from equation (2.5) after first-order integration, for m = 2, 3, i.e.
f (1) (2) i H (1) and f (2) By multiplying the first equation of (2.19) by D (1) while adding Ic 0 and IE inv(1) 1 leads to or alternatively as where D (2) = ( D (1) ) 2 and it can be proved that ξ = D (1) I. E inv(2) 2 is the error due to second-order iDQM-by-inversion operation. Repeating the same operation as in equations (2.21) and (2.22) twice for the second equation of (2.19), and adding the last three terms in equation (2.22), leads to (2.23) or alternatively as 3 and it can be proved that (1/2)ξ 2 = D (2) is the error due to third-order iDQM-by-inversion operation.
Equations (2.17) and (2.19) represent the second-and third-order equivalent of equation (2.8) from which a vector function F can be recovered from its higher derivatives. Accordingly, the general expression for the mth-order iDQM-by-inversion operation is given as If coordinate transformation is performed for the interval ξ ∈ [a b] to y ∈ [0 L] while noting from equation (1.3) that l i (0) = 0, equations (2.9) and (2.26) can be transformed to It should be noted that after transformation, point ξ p = y p = 0. Subsequent derivations will be presented in the new coordinates, y.
where l i (y), f (y i ) and E(f ) are described according to equations (1.2)-(1.5). Differentiating equation (2.29) once according to DQM approximation leads to the discretized equation, By evaluating the first derivative of the function y m analytically, equation (2.30) can be rewritten as which is further simplified (via inversion) after rearrangement in line with procedures in §d as which is now simplified as where c 0 is the intercept of f (y). Since f (y) = y m , such that the intercept c 0 = 0, then equation (2.35) becomes

Inverse differential quadrature method error
It is important to formulate an iDQM error estimate to describe the performance of iDQM in terms of numerical accuracy and numerical stability. Since iDQM-by-inversion is adopted in this work to circumvent computational issues arising from analytical or numerical integration, there is also a need to quantify the error estimate due to DQM inversion in order to assess the approximation quality of the proposed iDQM-by-inversion. Moreover, as the DQM coefficient matrix is used to compute iDQM weighting factors, a comparison of discrepancy between DQM and iDQM-byinversion error estimates is necessary to confirm improvements and drawbacks of the proposed method.
(a) Error formulation of inverse differential quadrature method-by-integration Consider a continuously differentiable function up to mth degree, f (y), where m is very high, Approximation of the first-order derivative of f , f (1) , according to equation (2.1) gives To determine the approximation error, a function, F(z), is defined such that So, by setting F(y) = 0, the approximation error is estimated as Applying Rolle's theorem [51] repeatedly leads to the Nth-order derivative of F(z), F N (z), which has at least one root, μ, between y 1 and y N leading to while noting that f In general, μ is a function of y.
In line with equation (3.7), let us denote the maximum error of f (1) approximation for a fixed order of y as E 1 max , which is expressed as where K 1 = |f N+1 (μ)| max . Now, the maximum error estimate of the original function, f (y), recovered via applying iDQM-by-integration to equation (3.7) is obtained thus (1) 1 max are the first-order integrals of M N (y) and E 1 max , respectively. Since K 1 and N are constants, it can be stated that the total error arising from numerical integration To have a clearer understanding of the implication of equation (3.9), consider a specific case of Chebyshev grid distribution, where y i are the roots of the Chebyshev polynomial, T k (y), of order k. In this context, y i can be expressed as (3.10) In connection with equation (3.10), the polynomial, M N (y), can be expressed in terms of T k (y) as (see [23]) is the first-order derivative of T N (y). Setting y = cosθ and T N (y) = cosNθ in equation (3.11) leads to, (3.12) The first-order integral of M N (y) in equation (3.12) is expressed as (3.13) Subject to equation (3.10), equation (3.13) can be recast in terms of y i at the grid points as: (3.14) As M int(1) N (y i ) varies in N by order O(N/N) in equation (3.14), it can be deduced that for a fixed order of y, the accuracy of iDQM operations is not substantially affected after first-order integration as N increases, since O(N/N) is mutually compensating (i.e. varies directly and inversely equally in order of N).
The error for recovering f from its second-order derivative, f (2) , can be obtained by redefining equation (3.3) as Then, repeating the procedures of equations (3.4)-(3.8) leads to the maximum error for f approximation from f (2) as Adopting a Chebyshev grid representation as applied in equations (3.10)-(3.14) results in A similar approach can be used to recover the maximum iDQM approximation error of f (y) from its third-order derivative, f (3) , as, (y) in the context of Chebyshev grid is given by (3.20) Equations (3.18) and (3.20) show that M . This observation illustrates that, for high-order approximation due to iDQMby-integration, the approximation error is scaled by a function which varies inversely in multiple orders of N. Therefore, subject to a fixed order of y, high-order iDQM-by-integration operation is potentially stable numerically as N increases.
(b) Error of inverse differential quadrature method-by-inversion According to equation (2.17), the approximation error for first-order iDQM-by-inversion reads Noting equation (3.9), the maximum approximation error for first-order iDQM-by-inversion becomes (0) to the total approximation error in equation (3.22) is expected to be minimal. In the same vein, the second-order approximation error according to iDQM-by-inversion is obtained according to equation (2.22), (1)  1 (0) in equation (3.24) using equations (3.9) and (3.17), the maximum approximation error due to second-order iDQM-by-inversion operation gives To establish a relationship between K 2 and K 1 , consider the Nth derivative of f (y) expressed as Furthermore, a recursive relation can be established between a derivative and its lower-order derivative as and Substituting equation (3.30) into equation (3.25), the approximation of the second-order iDQM inversion gives Applying similar procedures as in equations (3.24)-(3.32) leads to the approximation error for the third-order DQM inversion: According to equations (3.31) and (3.33), high-order approximation by iDQM-by-inversion incurs progressive error, which increases successively by order N.
(c) Comparison of inverse differential quadrature method error with differential quadrature method error As noted in Shu [23], the attributable error as a result of first-order DQM approximation is given by where μ is a function of y. Considering equation (3.8), the maximum error of function approximation by DQM is analogously expressed as where K 1 = |f N (μ)| max . By differentiating equation (3.35) to obtain the error due to first-order numerical differentiation, we get (1) max are the first derivatives of M N (y) and E max , respectively. According to equation (3.36), the error arising from first-order DQM differentiation is constrained by  (y), since K 1 and N are constants. Adopting Chebyshev grid representation, considering equations (3.11)-(3.12), the first derivative of M N (y) is expressed as Substituting for y i at the Chebyshev grid points while noting equation (3.10) results in In the same vein, the maximum error due to second-order DQM approximation is given as Substituting for y i in equation (3.40) at the Chebyshev grid points using equation (3.10) gives to the total errors subject to first-order numerical differentiation and numerical integration, respectively, affect the overall accuracy and numerical stability of the approximations. Comparing iDQM-by-integration and DQM error representation in this regard, it can be deduced that, for a fixed order of y, the resulting error from firstorder iDQM-by-integration is less than the first-order DQM approximation, as the total error magnification on account of the order of M dif(1) N and M int(1) N , respectively, in N is higher for DQM than iDQM approximation. This observation is also true for high-order iDQM-by-integration operation since, given a fixed order of y, the contribution of M int(1) N to the total error is less than the to the total high-order differentiation error of DQM operations. Although iDQM-by-inversion operation progressively magnifies the approximation error by order N, the rate of increase in the order of N for every successive inversion by iDQM operation is less than the rate of increase in the order of N for every successive differentiation by DQM. Therefore, iDQM operations in general potentially demonstrate superior numerical stability than DQM operations.

Numerical results and discussions
In this section, we present some illustrations of the numerical accuracy and numerical stability of the proposed iDQM for functional approximation as well as solution of ordinary differential equations (ODEs) and partial differential equations (PDEs) representing linear and nonlinear systems. We further demonstrate how different schemes of iDQM can be implemented using several examples. The results obtained are then benchmarked with DQM and exact solutions to evaluate the performance of iDQM solutions. To show the robustness of iDQM solutions, numerical analyses based on Lagrange basis polynomials on a non-uniform Chebyshev grid structure are performed and the errors computed. In addition, a comprehensive error analysis to examine the performance of iDQM schemes in terms of convergence and error propagation is performed using a fourth-order boundary value problem (BVP). The results are benchmarked against DQM and exact solutions to evaluate the gains of iDQM approximations.
Note: all examples performed in this work are implemented using iDQM-by-inversion since it proves more computationally efficient than iDQM-by-integration.
(a) Approximation of function and its higher derivatives According to table 1, iDQM estimates of f in equation (4.1) and its derivatives prove more accurate than DQM estimates. This observation is more evident for second-and third-order iDQM schemes in which the total error is less perturbed by numerical differentiation by DQM (to obtain high-order derivatives) or numerical integration by iDQM (to obtain low-order derivatives). Depending on the order of iDQM, the accuracy of functional approximations may fluctuate between low-order derivative approximation and low-order integral approximation, which is beneficial compared with DQM estimates, which accumulate error subject to successive highorder numerical differentiation. Apart from this, the rate of decline in numerical accuracy for DQM estimates is higher than for iDQM estimates, which is a good indication of numerical stability in favour of iDQM. This aspect is discussed further in §4g.
(b) Solution schemes for systems of differential equations by the inverse differential quadrature method By approximation of higher derivatives instead of the original function, the proposed iDQM formulation presents a unique opportunity to tune the order of a system, which aids in control of numerical accuracy and numerical stability of the system. In this context, two concepts are proposed in the following subsections. in a system, such that lower-order functions can be obtained via iDQM integration while higherorder functions are obtained via DQM differentiation. For example, a fourth-order system of equations can be represented by approximation of the second-order derivative in a MiDQM scheme, in which case the first-order derivative and original function are obtained by numerical integration while third-and fourth-order derivatives are obtained by numerical differentiation. This strategy leads to reduction in the DQM order required for the system solution and, by implication, reduction of the order of the DQM approximation error. This approach is highly promising, and therefore noteworthy, in that it allows tuning of the numerical accuracy of DQM to achieve improved solution. A demonstration of the implementation of this approach is presented for one dimension and two dimensions in appendix A.
(ii) Full inverse differential quadrature method In contradistinction to its DQM counterpart, FiDQM presents an opportunity to approximate the highest derivative in a system and then apply equation (2.28) to retrieve lower derivatives via iDQM operation. As established in the previous section, given a geometry, the error accrued by integrating a high-order function to get low-order estimates is quite stable numerically compared with differentiating low-order functions to get high-order estimates. As a result, depending on properties such as geometric specifications, boundary conditions or order of a system, high numerical accuracy can be achieved potentially by FiDQM schemes compared with DQM. A demonstration of the implementation of this approach is presented for one dimension and two dimensions in appendix A.  where E is Young's modulus of the beam and I is the second moment of area of the beam's crosssection. According to [52], the exact solution of the beam deflection is given by Given the iDQM discretization scheme described in appendix A, it is noted that the equations arising from iDQM constitute an underdetermined system because the number of unknowns exceed the number of equations. In this context, a pseudoinverse procedure based on truncated singular value decomposition described in [53] is adopted in this work to solve the systems of equations. According to figure 3, all iDQM estimates of the deflection and moment agree well with exact and DQM solutions showing the accuracy of iDQM solutions. In addition, according to table 2, the relative errors due to iDQM estimates computed for beam deflection, moment and shear force show significant improvement over DQM estimates for the same number of points. This remarkable improvement shows the great potential of iDQM schemes for numerical solution of ordinary differential equations.    Table 3. Maximum absolute error (relative error) of iDQM and DQM estimates for temperature and its derivatives.

(d) Nonlinear steady-state solution of heat conduction in slab with temperaturedependent conductivity (ODE)
The problem involves a steady-state heat conduction in a slab with temperature-dependent thermal conductivity in which the non-dimensional form of the temperature (θ) governing equation is expressed as The exact solution of the problem is given in [54] as The solution of the nonlinear system, i.e. equation (4.5), is obtained based on Newton-Raphson optimization after iDQM discretization. iDQM results prove accurate with respect to the exact solution and DQM estimates as temperature profile of the slab as well as its first-order derivative match satisfactorily ( figure 4). As expected, the error estimates reported in table 3 demonstrate that, although the accuracy of iDQM and DQM estimates of the temperature variable compares equally, iDQM estimates of higher-order derivatives of the temperature variable proves more accurate than DQM estimate, suggesting that error propagation during iDQM operations is less than for DQM operations.

(e) Solution of convection-diffusion equation (PDE)
Consider a steady-state convection-diffusion equation with boundary conditions as follows,   where P e is the Peclet number. The exact solution for the given partial differential equation is given in [55] as where σ = π 2 + P 2 e /4. According to figure 5, DQM and iDQM values converge to an exact solution of the PDE in equation (4.7) for a 41 × 41 grid, whereas in figure 5, only second-order iDQM converges to the exact solution for a 101 × 101 grid. It is quite evident from table 4 that the accuracy of the numerical values has a strong dependence on the Peclet number. For the case of Pe = 1000, DQM and first-order iDQM do not converge to the exact solution for the 101 × 101 grid and further grid refinement fails to improve the solution. According to the findings in [38], an upwind scheme is necessary to obtain accurate solutions for a high Peclet number. Nonetheless, second-order iDQM furnishes an accurate solution for Pe = 1000 without an upwind scheme, showing the accuracy of the proposed iDQM. Figure 6 demonstrates the agreement of second-order iDQM with an exact solution over the entire domain.  We consider a thin isotropic plate with dimensions a and b in x and y coordinates, respectively, simply supported on all the edges and under sinusoidally distributed load. The governing differential equation and associated boundary conditions are given as follows, where w is the transverse deflection of the mid-plane of the plate under loading q(x, y) = q 0 sin(π x/a) sin(π y/b), q 0 is the amplitude of the sinusoidally distributed load, D is the flexural stiffness of the plate given as D = (Eh 3 )/(12(1 − ν 2 )), E is Young's modulus, ν is Poisson's ratio and h is the thickness of the plate. Material and geometric properties of the plate are given by E = 200 GPa, ν = 0.3, a = b = 1 m, h = 0.01 m and q 0 = 1 Pa. The Navier's closed-form solution is simply given by, Solutions of the plate equation obtained by iDQM and DQM are shown in figure 7 where deflections and stresses for the plate match DQM and Navier's solutions demonstrating accuracy of iDQM. Figure 8 shows the deformed planform of the plate under the loading by fourth-order iDQM along with the maximum absolute error in the entire domain. Furthermore, the percentage error of DQM and iDQM estimates in table 5 clearly demonstrates faster convergence for fourth-, third-and second-order iDQM over DQM, highlighting the computational merits of the proposed method.

(g) Error analysis (measure of numerical accuracy)
To appropriately examine the convergence of iDQM, it is important to assess the numerical accuracy of iDQM estimates subject to increased discretization of the domain. In this regard, we (h) Error propagation (measure of numerical stability) As a measure of numerical stability, propagated error arising from approximation of low-order functions from high-order estimates (by using iDQM) or approximation of high-order functions from lower-order estimates (by using DQM) is examined in this section. Figure 10 shows error propagation of the BVP example in §4g, in which high-order functions (labelled superscript d) are computed from primary estimates (labelled superscript p) using DQM operation. On the other hand, low-order functions (labelled superscript i) are computed from primary estimates by iDQM-by-inversion. According to figure 10, the errors resulting from high-order integration by iDQM-byinversion are minimal, indicating numerical stability of iDQM operation. In the case of DQM, errors propagate by multiple orders for successive differentiation operation. While the error accumulation for low-order DQM approximations seems tolerable, the multiple order increment in the error due to high-order differentiation can cause inaccuracy which, in turn, leads to  numerical instability. As already mentioned in §3c, numerical stability of the DQM solution is significantly affected by the high-order derivative of M N (y), which increases geometrically as the order of N increases. Therefore, on increasing N, high-order derivatives of M N (y) quickly offset the total error to reach a lower error bound. However, the lower error bound of low-order estimates from high-order estimates remains stable after iDQM-by-inversion, indicating low error propagation and improved numerical stability.

Computational efficiency of the inverse differential quadrature method
To measure the computational efficiency of iDQM, we consider a comparison of the time and space memory requirements for the different benchmark problems to converge to a fixed value of the maximum absolute error max for DQM and iDQM estimates. In this regard, the bandwidth (b) and the total primary degrees of freedom (n) of the final matrix A for a given algebraic system Ax=b are computed according to table 6.
Some benchmark examples (4c, 4d and 4g) are chosen in table 6 to reflect different types of analysis, boundary conditions and nonlinearities that directly affect the computational complexities of a given numerical problem. According to table 6, iDQM approximation preserves the order of the numerical complexities of the problem in terms of time and space requirements as DQM. It is worth noting that the DQM estimate for example (4c) fails to converge to the threshold absolute maximum error max , i.e. 10 −13 . Thus, it is concluded that iDQM approximation preserves the computational efficiency of DQM approximation.

Conclusion
This study proposes a novel iDQM for numerical analysis of engineering systems. Given a system of high-order differential equations, the proposed iDQM approximates high-order variables rather than the original function which can be subsequently recovered by integration. To deal with issues bordering on computational inefficiency of analytical integration and numerical complexity of Gaussian integration, this study develops a novel strategy which relies on inversion of the existing formula of the DQM to compute the required iDQM weighting factors. Furthermore, to evaluate the performance of iDQM solutions, detailed derivations of iDQM error estimates based on integration and DQM inversion are developed in this work, which are then compared with DQM error estimates outlined in [23]. In the context of the iDQM scheme, two implementation approaches identified as Mixed iDQM and Full iDQM are proposed to obtain solutions of the examples provided in this work. Remarkably, the concept of Mixed iDQM provides an excellent opportunity to control the accuracy of system solutions by combining the numerical advantages of low-order differentiation and low-order integration to achieve an improved solution. Subsequently, a demonstration of iDQM implementation for functional approximation, and numerical solutions of systems of high-order ordinary differential equations and partial differential equations representing linear and nonlinear systems, prove that iDQM operations are potentially robust to furnish accurate solutions to numerical systems. Finally, an appraisal of the convergence and numerical stability of the iDQM approach suggests that, compared with DQM, improved convergence can be obtained for systems solution and improved numerical stability is guaranteed by using the proposed method without loss of computational efficiency. is a contributing author who also supervised the project, and contributed to numerical analysis. All authors gave final approval for publication and agree to be held accountable for the work performed therein.