A bivariate Chebyshev polynomials method for nonlinear dynamic systems with interval uncertainties

A bivariate Chebyshev polynomials approach is proposed to estimate the dynamic response bounds of nonlinear systems with interval uncertainties. The existing collocation method directly searches the maximum and minimum values of the surrogate model in the entire interval space by the scanning method (SM). The presence of too many uncertain parameters will lead to expansive computational cost. To overcome this shortcoming, the dynamic response is decomposed by a bivariate function decomposition (BFD), established based on high-order Taylor expansion, into the sum of multiple univariate and bivariate response functions. The above univariate and bivariate functions are fitted using Chebyshev polynomials, and polynomial coefficients are obtained through one-dimensional (1D) and two-dimensional (2D) interpolation points. Thus, the solution of the nonlinear dynamic systems with uncertain parameters can be transformed into that of univariate and bivariate Chebyshev interval functions. The extremum values of the low-dimensional Chebyshev interval functions can be found by SM, and then the bounds of dynamic response are acquired by interval arithmetic. Since SM searches for extreme values only in 1D and 2D uncertain domains, the amount of calculation is reduced compared to searching the whole uncertain space. The efficiency, practicability and effectiveness of the proposed interval uncertainty analysis method are proved by three dynamic examples.


Introduction
Dynamic response forecasting of nonlinear systems is of great significance for structural design and analysis of most practical engineering. Traditional dynamic analysis is based on the deterministic model. However, any dynamic system has inevitable uncertainties, such as load fluctuation, manufacturing instability and dispersion of material properties, which make structural parameters uncertain [1][2][3][4]. Hence, the dynamic response of nonlinear systems is also uncertain. Noticeably, small changes associated with parameter uncertainty may largely change the dynamic response [5]. As a result, the influence of uncertainty on dynamic systems must be considered.
In uncertainty analysis, probabilistic and nonprobabilistic methods are often employed to quantify the effect of uncertainty on dynamic response. In general, the probability method is the most valuable solution under the premise that the probability density function (PDF) of uncertain parameters is clearly defined [6][7][8]. Unfortunately, the method has an obvious disadvantage, as it needs complete statistical data to accurately describe the PDF of random variables or fields, and these data are usually inaccessible and expensive in engineering problems. Under the circumstances, some non-probabilistic methods, including fuzzy set [9], evidence theory [10] and interval set [11], are developed to describe the uncertainty caused by limited data and knowledge.
In interval sets, uncertainty is represented by the upper and lower limits of uncertain parameters, rather than other stochastic data. Therefore, it is an effective assistance to the conventional probability method. In the past two decades, interval sets and interval analysis approaches have been widely used in structural dynamic response investigation. The interval arithmetic (IA), which was the first applied for dynamic response analysis, directly obtains the response interval based on interval algebraic operation rules [12][13][14]. Since the correlation between uncertain parameters is neglected, IA usually leads to overestimation of interval values, which is called the wrapping effect. This effect is more evident and serious in multidimensional and non-monotonic dynamic systems. There is extensive research on this topic [15][16][17]. Besides, the affine algorithm [18][19][20] introduced to remedy the wrapping effect can alleviate overestimation, but can cause serious errors in the case of highly nonlinearity. This is because nonlinear affine operation will induce more noise terms and overestimate response bounds [20]. As long as the uncertainty problem is monotone, the vertex method can get accurate interval results [21][22][23]. However, it can hardly guarantee or prove monotonicity, especially in the field of structural dynamics.
The Taylor expansion-based method (TEBM) performs the first-or second-order Taylor series expansions at the midpoint of uncertain parameters and obtains the dynamic response interval [2,[24][25][26]. The advantage of TEBM is high computational efficiency, and the calculation consumption is raised linearly as the number of interval parameters rises. Meanwhile, it can reduce the wrapping effect due to the dependence among multiple interval parameters. However, TEBM has some disadvantages. First, it is accurate only when the range of uncertain parameters is narrow enough. The subinterval approach is a solution to this problem [27][28][29]. It divides the large interval into several small intervals, which improves the accuracy of interval analysis. But the amount of calculation will sharply increase as the number of subintervals increases. Second, TEBM requires the first-or even higher-order partial derivatives of the dynamic response concerning the interval variables. In some engineering conditions, the explicit functional relationship between input and output is too complicated, which may make it impossible to acquire first-and higher-order derivatives. Most importantly, TEBM can only moderately relieve overestimation, but may even lead to worse results in case of strong nonlinearity [25,30].
The dynamic response is always highly nonlinear relative to interval parameters. Compared with the above methods, the collocation method (CM) can solve large uncertainty and strongly nonlinear problems. Firstly, CM establishes a surrogate model according to the selected interpolation points. The proxy models mainly include Chebyshev polynomials [31][32][33], Bernstein polynomials [34], Legendre polynomials [35], Fourier series [36][37][38], arbitrary polynomial chaos [11], Szász-Mirakyan operator [39,40] and Baskakov operator [41]. Among these surrogate models, Chebyshev polynomials have been widely used in interval analysis of dynamic systems. Then, the auxiliary algorithms are applied to find the extreme values of the surrogate model, and obtain the interval results of the dynamic response. The auxiliary algorithms common include IA [32,42], genetic algorithm (GA) [43][44][45] and SM [35,[46][47][48]. IA has the fastest computing velocity, but often yields calculation errors due to the wrapping effect, especially when the dynamic response is largely uncertain and nonlinear. To this end, many improved IA-based methods have been studied, such as trigonometric function [49] and modified affine basis [50]. These methods only moderately improve the accuracy, but large errors still exist with the precise dynamic interval. Generally, GA achieves the predetermined accuracy through iteration, and its computation is expensive if there are many uncertain parameters. Moreover, when the nonlinearity of the surrogate model is strong, GA may fall into a local optimum. Compared with IA and GA, SM is the most precise, but its samples rise exponentially with the increased number of interval parameters.
This study concentrates on the following issues: When SM is employed to compute the maximum and minimum values of the surrogate model, the existing collocation methods need to search the entire uncertain domain. Thus, when there are many interval parameters, it easily falls into the curse of dimensionality. To increase efficiency and ensure accuracy, this paper proposes a bivariate Chebyshev polynomials method for interval analysis of dynamic response with many uncertain parameters. In this method, BFD decomposes the original dynamic function into the sum of univariate and bivariate response functions, which are then fitted by Chebyshev polynomials. In this way, the nonlinear dynamic systems described by ordinary differential equations (ODEs) can be transformed into 1D and 2D ODEs, and the whole interval space is divided into a linear combination of 1D and 2D spaces. Moreover, the ODEs can be solved by an appropriate numerical solver. Finally, the upper and lower bounds of the dynamic response at each time step are determined using SM and IA.
The statement of nonlinear dynamic system with uncertain parameters is given in Sect. 2. Section 3 illustrates the detailed derivation of Chebyshev interval method based on IA or SM. Then, a bivariate Chebyshev polynomials method is proposed in Sect. 4. In Sect. 5, the three dynamic cases are presented to verify the effectiveness and practicability of the proposed method.

Overview of nonlinear dynamic system with interval parameters
In general, the nonlinear dynamic system is represented by a group of ODEs. In the case of l state vectors, the l-dimensional ODEs are demonstrated as: where g ¼ g 1 ; g 2 ; . . .; g l ½ T is a vector containing nonlinear functions, y ¼ y 1 ; y 2 ; . . .; y l ½ T is a dynamic response state vector and y 0 is the initial condition vector; x ¼ x 1 ; x 2 ; . . .; x n ½ T is an n-dimensional vector; t t 0 ; t e ½ is the time range. Equation (1) is a deterministic model that can be computed by any numerical integration approach to obtain the change of y over time.
In actual engineering, the physical and geometric properties of the structure are uncertain to varying degrees. For uncertain input vector x, if the upper and lower bounds of the parameters can be estimated, then x can be expressed as [24]: where the superscript I is the interval; x C k and Dx k represent the midpoint and radius of the kth uncertain parameter, respectively; the transformation parameter x C k and Dx k are defined as: where x U k and x L k represent the upper and lower bounds of x k , respectively. Since the parameter uncertainty is transitive, the original ODEs become a group of interval ODEs. Therefore, the solution of Eq. (1) subjected to Eq. (2) is also an interval set [3]: Generally, the set C has a complex region. In interval mathematics [12], solving the nonlinear problem of Eq. (1) constrained by Eq. (2) is synonymous to find a multidimensional rectangle defined by Eq. (4). In other words, the upper and lower bounds of Eq. (4) will be calculated as follows: and Obviously, the interval value of the dynamic response is essentially the supremum and infimum of Eqs. (6) and (7). The simplest and most direct method to solve Eq. (5) is SM [32], but it is computationally intensive. Additionally, IA and TEBM are simple, effective and commonly applied to evaluate interval functions. However, they may lead to overestimation as mentioned in the introduction. Moreover, TEBM is only suitable for small uncertainty problems. To solve the dynamic response problem with large uncertainty and strong nonlinearity, the Chebyshev interval analysis method is proposed, as elaborated in Sect. 3. Because of the expensive computation in scanning the whole interval space, the proposed method will solve this problem accurately and effectively in Sect. 4.

Improved Chebyshev interval model
With a combination of Chebyshev polynomials of degree p, the dynamic response y x ð Þ in the uncertain domain X ¼ À1; 1 ½ n can be approximately expressed as [31,32]: where h is the number of zeros appearing in the subscript i 1 . . .i n ; b i 1 ...i n is the expansion coefficient; . . .; x n and X ¼ ½X 1 ; X 2 ; . . .; X n are transformed as follows: Let T p be the Chebyshev polynomial of degree p, then [51]: For n-dimensional problems, T i 1 ...i n X ð Þ is the tensor product of each one-dimensional Chebyshev polynomial T p X i ð Þ: In addition, each Chebyshev polynomial T i k X k ð Þ satisfies the following orthogonal relationship: where iis the inner product of the weight function in the specified domain ½À1; 1 n . Equation (12) can be further expanded: According to the orthogonality of Eq.
Given the difficulty to directly solve the integral formula of Eq. (14), the Gauss integration technique [52] where n j k is the interpolation point. The zeros of the Chebyshev polynomial of degree (p ? 1) are the interpolation points, which can be computed as: The number of interpolation points M for an ndimensional problem is M¼ p þ 1 ð Þ n from Eq. (16). Obviously, M will be large for high-dimensional problems, leading to expensive calculation. To save costs, only a part of the interpolation points are utilized to build the Chebyshev polynomial. Therefore, the Chebyshev surrogate model in Eq. (8) is improved to y IC x ð Þ as [31]: In Eq. (17), the terms above p are deleted, and k ¼ n þ p ð Þ!= n!p! ð Þ terms are left. Since not all interpolation points are applied to construct y IC x ð Þ, the coefficients in Eq. (15) cannot be calculated. Thus, the least squares method is used to generate the coefficients, and the rules for selecting interpolation points can be described as [53]: Considering the interval uncertainty in Eq. (17), the improved Chebyshev interval model (IC) of nonlinear dynamic response can be expressed as: After that, the bounds of y IC x I ð Þ are computed. Two methods are commonly used: (1) improved Chebyshev interval model using interval arithmetic (IC ? IA); (2) improved Chebyshev interval model using scanning method (IC ? SM).

Improved Chebyshev interval model using interval arithmetic
Since can be expressed as: By expanding Eq. (20), the lower bounds of y IC IA x I ð Þ can be stated as: The upper bounds of y IC IA x I ð Þ can be defined as: Equations (21) and (22) can reduce the overestimation to a certain extent, but cannot completely eliminate the overestimation [32].

Improved Chebyshev interval model using scanning method
To reduce the error introduced in solving the bounds of cos i 1 h I 1 Á Á Á cos i n h I n using IA, the interval value of First, all the sampling points enclosed in the interval x I are selected. In each interval the set of all sampling pointsx is formulated as: In Eq. (23), there are a total of e n sampling points. Second, each given sampling point is substituted into Eq. (17) to calculate the value of y ICx ð Þ. All values of y ICx ð Þ form the set are: Third, by comparing the values of the elements in C IC , the lower bounds of y IC x I ð Þ can be obtained: Similarly, the upper bounds of y IC x I ð Þ is: Equations (25) and (26) are applied to scan the explicit IC constructed by Eq. (19), which improves the computational efficiency compared with the direct scanning of the original dynamic response function. However, considering that the sampling points of IC ? SM is e n , the calculation cost is still unacceptable for the dynamic response problem with more sampling locations or more interval parameters.

Construction of Bivariate Chebyshev interval model
To reduce the calculating consumption of scanning the IC, we first decompose the n-dimensional dynamic response function into a sum of univariate and bivariate response functions by BFD, and then fit these functions based on Chebyshev polynomials.

Bivariate function decomposition
The original dynamic response function y x ð Þ is processed by high-order Taylor expansion at the midpoint x C of the uncertain parameters, and then retained to the second-order term [54]: Similarly, these univariate points If we consider (27), so y x ð Þ becomes a function containing x 1 and x 2 : Similarly, these bivariate points are considered in Eq. (27) separately: Second, adding each term in Eqs. (28) and (29), the sum of univariate function is expressed as: Similarly, the sum of the bivariate function can be obtained by adding Eqs. (30) and (31): By shifting the terms of Eqs. (32) and (33) separately, we can get: and X 1 Finally, substituting Eqs. (34) and (35) into Eq. (27), the bivariate function decomposition (BFD) of y x ð Þ can be expressed as: If y y C 1 ; . . .; y k 1 ; . . .; y k 2 ; . . .; y C n À Á , y x C 1 ; . . .; x k ; . . .; x C n À Á and y x C ð Þ are abbreviated as y x k 1 ; x k 2 ð Þ, y x k ð Þ and y 0 , respectively. Thus, Eq. (36) is rewritten as: where y x ð Þ is comprised of n univariate functions and n n À 1 ð Þ=2 bivariate functions. Each univariate function and bivariate functions contain only one and two uncertain parameters, respectively, and the other parameters are replaced by their midpoint values.

Bivariate Chebyshev interval model
The above y x k ð Þ and y x k 1 ; x k 2 ð Þ are fitted by Eq. (8) separately: and where b i k and b i k 1 i k 2 represent the expansion coefficient with one uncertain parameter x k and two uncertain parameters (x k 1 and x k 2 ), respectively. According to the judgment rules given by Eq. (18), each univariate and bivariate functions should select (p ? 1) and (p ? 1) 2 interpolation points, respectively. Thus, b i k and b i k 1 i k 2 can be solved by Eq. (15). Substituting Eqs. (38) and (39) into Eq. (37) and considering interval uncertainty, the bivariate Chebyshev interval model (BC) can be expressed as: where y x I k À Á and y x I k 1 ; x I k 2 are the univariate and bivariate Chebyshev interval function, respectively, which can be calculated by: The interpolation points required by the Chebyshev model (CM, Eq. (8)), IC, and BC are given in Table 1.
The number of interpolation points of CM increases exponentially with the rising number of uncertain parameters, but those of IC and BC are dramatically reduced, avoiding the huge calculation costs (Table 1). Furthermore, BC has fewer interpolation points for high-dimensional problems.
After BC is obtained as shown in Eq. (40) denote a 1D sampling point, then the set of all 1D sampling pointsx k is: Thus, each univariate function requires e sampling points.
For y x k 1 ; x k 2 ð Þ, e equidistant sampling positions are also selected in the intervals x I k 1 and x I k 2 . Let sampling point, then the set of all 2D sampling points is:^x Therefore, each bivariate function requires e 2 sampling points. The set of all sampling pointsx BC can be stated as: Obviously, scanning BC requires n(n-1)e 2 /2 ? en sample points in total.
Second,x k is substituted into Eq. (38) Third, by comparing the values of the elements in C k and C k 1 k 2 , the lower bounds of y x I k À Á and y x I can be expressed as: Similarly, the upper bounds of y x I k À Á and y x I k 1 ; x I k 2 can be described as: Finally, substituting Eqs. (47) and (48) into Eq. (41), the lower bounds of BC can be obtained based on IA as: Similarly, the upper bounds of BC can be expressed as: Compared with IC ? SM in Eqs. (25) and (26), which requires e n samples, BCM only needs n(n-1)e 2 / 2 ? en. The sample distribution comparison between IC ? SM and BCM in a three-dimensional uncertain space is shown in Fig. 1. Clearly, BCM has a relatively low computational cost (Fig. 1).

BCM to solve nonlinear dynamic systems with interval parameters
BCM is introduced for nonlinear dynamic systems expressed by ODE with interval uncertainty. When employing numerical methods to solve differential equations, the time range t t 0 ; t f Â Ã is discretized into a group of time points t 0 ; t 1 ; . . .; t b ; . . .; t f À Á , and the solution at each time point is computed. At time t b , the dynamic response concerning interval parameters can be approximated by BC: In Eq. (51), the original dynamic response is approximated as the sum of multiple univariate and bivariate functions. Then, these functions are explicitly fitted by Chebyshev polynomials, and the expansion coefficients are obtained through solving 1D and 2D ODEs. When the surrogate model is built, the interval of y at time t b can be determined simply by auxiliary methods. Moreover, any numerical method can be used to solve 1D and 2D ODEs. The flowchart of the solution is depicted in Fig. 2: Step (1). Define the degree p of the Chebyshev polynomial, the interval of x, the sampling position e of each interval parameter, the length of the time period t f , and the time step Dt.
Step (2). Decompose the original dynamic response into the sum of multiple univariate and bivariate functions by BDE (Eq. (36)); Step (3). Generate 1D and 2D interpolation points by Eq. (16) for fitting the above univariate and Step (4). Substitute the 1D and 2D interpolation points separately into the original differential equation, and use the numerical method to compute the actual 1D and 2D dynamic responses at the time t b : y n j k ; t b À Á and y n j k 1 ; n j k 2 ; t b ; Step (5). Calculate the 1D and 2D expansion coefficients, b i k and b i k 1 i k 2 , to build univariate and bivariate Chebyshev interval functions at the time t b ; Step (6). Insert the 1D and 2D sampling points into the univariate and bivariate Chebyshev interval functions, and obtain the sets C k and C k 1 k 2 ; Step (7). Scan to find the maximum and minimum values of C k and C k 1 k 2 , and get the bounds of univariate and bivariate Chebyshev interval functions; Step (8). Use IA to compute the interval value of y at the time t b ; Step (9). If b [ t f Dt, output the interval change of y over time; otherwise, let b ¼ b þ 1, and return to step 4.

Numerical cases
Three numerical cases are used to verify the effectiveness and accuracy of the proposed method: (1) a double pendulum problem with three interval parameters, (2) a damped spring-mass system with five interval parameters, and (3) a two-degree-of-freedom (DOF) vehicle model with seven interval parameters. SM (regarded as the exact solution), IC ? IA, IC ? SM and BCM are used for comparison. All the algorithms are programmed using MATLAB R2019a operated on the same computer. The numerical calculation adopts the Runge-Kutta method.

Double pendulum problem
In this double pendulum system (Fig. 3), h 1 and h 2 are the angles of the pendulum, m 1 and m 2 are the masses of the two rods, L 1 and L 2 are the lengths of the two rods, and g is the local acceleration of gravity.
The ODEs of the highly nonlinear double pendulum system are expressed as [22]: where w 1 and w 2 represent the angular velocity of the two pendulum rods, respectively. In this example, L 1 , .81] m/s 2 , respectively. The initial condition is Fig. 2 the solution process of the BCM for nonlinear dynamic systems with interval parameters The SM, IC ? IA, IC ? SM and BCM are used to solve ODEs with interval parameters at time of 0-8 s. The time step is Dt ¼ 0:001 s. Moreover, the degree p of Chebyshev polynomial is 3, and each interval parameter uses 20 symmetrical sampling positions. The calculation results are plotted in Figs. 4 and 5. For intuitive comparison, ten time points are chosen, and the details are listed in Tables 2 and 3. IC ? SM and BCM provide more accurate results than IC ? IA, and IC ? SM and BCM almost consistent with SM ( Figs. 4 and 5, Tables 2 and 3). Moreover, the results of IC ? IA diverge slowly with time. In contrast, IC ? SM and BCM can effectively reduce the wrapping effect throughout the numerical process.
Compared with SM, IC ? SM and BCM are both searching for the maximum and minimum values of the surrogate models, and thus take less computational consumption. IC ? SM requires 20 3 = 8000 sampling points at each time point for the double pendulum problem with three interval parameters. However, BCM only needs 1261 sample points, which is more efficient than IC ? SM.

Damped spring-mass system
In this two-DOF damped spring-mass model (Fig. 6), M 1 and M 2 are masses, C is the damping coefficient,  The ODEs of the system can be expressed as [55]: where y 1 and y 2 are displacements of M 1 and M 2 , respectively. The initial condition is _ y 1 ; _ y 2 ; y 1 ; y 2 ½ t¼0 ¼ 0; 0; 0:1; 0:1 ½ j . Here, M 1 , M 2 , C, K 1 and K 2 are regarded as uncertain parameters (Table 4). SM, IC ? IA, IC ? SM and BCM are employed to solve ODEs with interval parameters in 0-100 s. The time step is Dt ¼ 0:02 s. Given the degree p of the Chebyshev polynomial is 2, each interval parameter uses ten symmetrical sampling positions. Figures 7  and 8 show the interval bounds of y 1 and y 2 in the time period 100 s, respectively. Moreover, ten time points are selected, and the details are listed in Tables 5 and 6.  Fig. 7 The interval bounds of the displacement y 1 obtained by four interval methods

Fig. 6 A damped spring-mass system
The y 1 obtained by IC ? IA, IC ? SM and BCM has a smaller error in the linear region (i.e., in the red dashed box) (Fig. 7), which can also be illustrated by y 2 in the red dashed box of Fig. 8. However, in the highly nonlinear region (the peaks and valleys of the curve), the error between the dynamic response intervals obtained by IC ? IA and SM is considerably large. The displacement intervals of IC ? SM and BCM are always close to the accurate solution without major deviation (Figs. 7 and 8, Tables 3 and 4). Therefore, IC ? SM and BCM can acquire correct and dependable dynamic response interval throughout the numerical period.
For the damped spring-mass system with five interval parameters, IC ? SM requires 10 5 sample points at each time point, while BCM only needs 1051 sample points, so the computational efficiency is much higher than IC ? SM.

Two-DOF vehicle model
To prove that the proposed method is applicable to practical engineering, the third example adopts a two-DOF vehicle model ignoring the influence of steering and suspension system (Fig. 9). In Fig. 9, O is the center of mass; v and u are the lateral and longitudinal velocities of the center of mass, respectively; b is the sideslip angle of the center of mass; w r is the yaw velocity; a and b are the distances from the center of mass to the front and rear axles, respectively; d is the front wheel steering angle.
The ODEs of the two-DOF vehicle can be written as [56]: where k 1 and k 2 are the front and rear wheels cornering stiffness, respectively; m is the vehicle mass, I z is the moment of inertia. b can be denoted as:  Fig. 8 The interval bounds of the displacement y 2 obtained by four interval methods Substituting Eq. (55) into Eq. (54), we can get: where This example considers the dynamic response of the steering angle step input: The step input is defined as the front tire rotating to 2°in 0.2 s, and then keeping the angle unchanged, which can be expressed as: d¼ 10p 180 t; t\0:2 Assuming that all parameters of the vehicle model are uncertain, the information of uncertain parameters is listed in Table 7. The midpoint values of these parameters are all from Ref. [56]. Moreover, the initial condition is _ w r ; w r ½ t¼0 ¼ 0; 0 ½ j .  Clearly, this dynamic model is relatively complex. SM, IC ? IA, IC ? SM and BCM are all applied to solve this dynamic systems in a time period of 0-3.5 s. The time step is Dt ¼ 0:001 s. Given the degree p of the Chebyshev polynomial is 2, each interval parameter uses ten symmetrical sampling positions. The interval bounds of w r are shown in Fig. 10. Besides, ten time points are chosen ( Table 8).
The interval bounds of w r generated by IC ? SM and BCM are basically consistent with that of SM, while IC ? IA has a large error (Fig. 10). Table 8 and the red dotted box in Fig. 10 demonstrate this result more intuitively. Therefore, BCM has a good accuracy even for complex dynamic systems. However, due to the use of IA in BCM, some overestimation cannot be completely avoided.
Given that IC ? SM needs 10 7 sample points at each time point for the vehicle model with seven interval parameters and that BCM only requires 2171 samples, BCM has a very high computational efficiency. Furthermore, the example clearly shows that BCM can achieve high accuracy and low calculation amount when computing the response limits of a nonlinear dynamic system with interval uncertainty.

Conclusions
A bivariate Chebyshev polynomials method is put forward to estimate the upper and lower bounds of nonlinear dynamic systems with interval parameters. This method can effectively alleviate the curse of dimensionality caused by scanning the entire uncertain domain. Based on BFD and Chebyshev polynomials, the solving of high-dimensional ODEs with interval parameters is converted into that of 1D and 2D ODEs, and the entire uncertain space is decomposed into a linear combination of 1D and 2D spaces. The upper and lower bounds of the dynamic response at each time step can be easily approximated by SM and IA. Three numerical examples are used to illustrate the feasibility and effectiveness of the method, and the following conclusions can be further summarized: BCM provide more precise and reliable dynamic response limits than IC ? IA. Meantime, the interval value of BCM is almost consistent with those of IC ? SM and SM in the whole numerical period, but its computational efficiency is much higher. Hence, BCM is more accurate with less computational cost in solving nonlinear dynamic systems with interval parameters. BCM is also non-intrusive and only requires 1D and 2D interpolation points to obtain a surrogate model of the original response. Thus, BCM is more easy-to-use in practical engineering, particularly for implicit dynamic response with multiple uncertain parameters.