Milling Stability Prediction: A New Approach Based on a Composited Newton–Cotes Formula

Based on a composited Newton–Cotes formula, this paper proposes a numerical method to predict milling stability considering regenerative chatter and focusing on rate and prediction accuracy. First, the dynamic model of milling motion is expressed as state-space equations considering regenerative chatter, with the tooth passing period divided into a set of time intervals. Second, a composited Newton–Cotes formula is introduced to calculate the transition function map for each time interval. Third, the state transition matrix is constructed based on the above-mentioned transition function, and the prediction stability boundary is determined by the Floquet theory. Finally, simulation analysis and experimental verification are conducted to verify the effectiveness of the proposed method. The simulation results demonstrate that, for the milling model with a single degree of freedom (DOF), the convergence rate and prediction accuracy of the proposed method are higher than those of the comparison method. The experimental results demonstrate that, for the milling model with two DOFs, the machining parameters below the prediction stability boundary can avoid the chatter as much as possible, ensuring the machined surface quality.


Introduction
Regenerative chatter arising from a self-excitation chatter vibration may be unavoidable in milling processes, which will reduce machining quality and efficiency and accelerate tool wear [1,2]. Researchers have attempted to avoid regenerative chatter using prediction [3,4], identification [5][6][7], and suppression [8][9][10] methods. The research results demonstrate that selecting machining parameters to avoid chatter occurrence is vital and that the stability lobe diagrams provide a reference for appropriately selecting machining parameters [11]. The objective of this study is to develop a method to draw stability lobe diagrams with high precision.
As regenerative chatter is integrated into the dynamics model, the dynamics model is reformulated as delay differential equations (DDEs) with time-periodic coefficients [12,13]. So far, many methods, including experimental, analytical, and numerical methods, have been proposed to solve the DDEs to obtain stability lobe diagrams. Based on the zerothorder harmonics of the Fourier series expansion, Altintas et al. [14] transformed the DDEs into the Laplace frequency, which is more suitable for high-radical-immersion milling. Taking higher-order harmonics into account, Merdol et al. [15] proposed a multi-frequency solution for low-radical-immersion milling, which is also called the multi-frequency method. In addition to the above-mentioned stability analysis methods in the frequency domain, there are some stability analysis methods in the time domain. Insperger et al. [16,17] first presented a semi-discretization method (SDM) for a delayed system and then presented an updated version of the SDM (first-order SDM) for periodic systems with a single discrete time delay. By extending the SDM and first-order SDM, the second-order SDM (second-order SDM) has also been proposed for variable spindle speed milling [18]. Ding et al. [19] proposed a full-discretization method (FDM) in virtue of a direct integration scheme. To further improve the prediction accuracy, Ding et al. [20] and Yan et al. [21] presented the second-order FDM and third-order FDM, respectively. Based on the integral equation and numerical integration formulas, Ding et al. [22] developed a concise semi-analytical method for constructing a discrete dynamical map approximating the DDEs. Dong et al. [23] have recently proposed an updated numerical integration method. Li et al. [24] proposed an accurate and fast milling stability prediction approach based on the Newton-Cotes method (NCM). In addition to the above-mentioned methods, other methods such as the Adams-Moulton-based method [25], the spectral element method [26], and the Lyapunov-Krasovskii method [27] exist. The aforementioned stability analysis methods mainly focus on three targets: convergence rate, prediction accuracy, and computational efficiency. These three targets usually contradict each other to a certain extent. Based on a composited Newton-Cotes formula, this paper presents a novel milling stability prediction method, i.e., the composited Newton-Cotes method (CNCM), with a high convergence rate and prediction accuracy.
The remainder of this paper is organized as follows. Section 2 provides a brief description of the composited Newton-Cotes formula. The composited Newton-Cotes-based milling stability analysis method is presented in Section 3. Comparisons of the convergence rate, prediction accuracy, and computational efficiency between the proposed method and NCM are presented in Section 4, and experimental verification is presented in Section 5. Finally, conclusions are drawn in Section 6.

Composited Newton-Cotes Formula
The nth-order Newton-Cotes formula has at least n-order algebraic accuracy [28], which is defined as where [a, b] denotes the parameter interval, and it is divided into n equal intervals. The span of each interval is defined as (b − a)/n, and then x k = a + k(b − a)/n is deduced. The C (n) k denotes the cotes coefficient and is defined as where depends only on n. It is worth mentioning that the approximation accuracy of Equation (1) is significantly reduced when n is greater than or equal to eight. Increasing n cannot ensure the approximation accuracy of Equation (1). To this end, the composited Newton-Cotes formula is introduced below. To improve the integration accuracy, midpoints of the interval should be considered. So the interval [x k , x k+1 ] is divided into four equal sub-intervals: x k , x k+1/4 , x k+1/2 , x k+3/4 , and x k+1 . The lower-order Newton-Cotes formula is used in each sub-interval, and then the approximation results of each subinterval are finally added together. According to Bool's rule (n = 4), the fourth-order composited Newton-Cotes formula [29] can be introduced as follows: where ∆ = (b − a)/n. The remainder of the equation above is

Mathematical Model
In this work, it is assumed that the workpiece is rigid, that the milling cutters have two DOFs and that there is no tool run-out situation in the milling process. The tool helix angle is not considered either. Considering the regenerative chatter, the mathematical model of milling motion can be formulated as a second-order differential Equation (24) as follows: where M, C, and K denote the modal stiffness, mass, and damping matrices of the milling system, respectively; q(t) denotes the displacement vector, . q(t) and .. q(t) denote its firstand second-order derivatives with respect to t; a p denotes the depth of cut; K C (t) is a time-periodic coefficient matrix that satisfies K C (t) = K C (t + T); T denotes the time delay. For the milling cutter, it is the same as one-tooth passing period, i.e., T = 60/(NΩ), where N denotes the number of tool teeth and Ω denotes the spindle speed (rpm).
A new vector x(t) as x(t) = q(t), .

q(t)
T is defined, and then the state-space expression of (5), can be expressed as [22] . where where K t and K n denote the tangential and normal linearised cutting force coefficients, respectively. The φ j (t) denotes the angular position of the jth tooth. Based on the Volterra integral equation of the second type [30], the analytical expression of (6) is deduced as follows: where t st is the starting time of the computation. Depending on whether the cutter touches the workpiece, the vibrations can be classified into free and excited vibration process [31]. If the cutter is not in contact with the workpiece, the vibrations are free vibration processes, which are equivalent to B(t) = 0. Then, Equation (11), can be simplified as where t f is the free vibration cycle. If the cutter is in contact with the workpiece, the vibrations are excited vibration processes. The corresponding cycle t i is where i = 2m m = 0, 1, 2, · · · , ceil(n/2) − 1, where ceil (·) is a function that rounds positive numbers to plus infinity) and h = T − t f /n. By substituting Equation (13), into Equation (11), we obtain Based on the composite Newton-Cotes formula introduced in Section 2, if the time interval [t i , t i+1 ] is divided into four equal sub-intervals, that is, t i , t i+1/4 , t i+1/2 , t i+3/4 , and t i+1 , then the state items x(t i+1/4 ), x(t i+1/2 ), and x(t i+3/4 ) can be approximated as follows: To simplify the process, several abbreviated expressions are used as below. Let By combining Equations (15)- (20), x(t i+1 ) can be computed by On this basis, we utilise Simpson's rule to calculate x i+2 : By combining Equations (21) and (22), we can obtain general iterative algorithms, and the state transition equations can be written as where where I is the identity matrix. After solving the equations above, the state transition matrix is established as follows: Micromachines 2023, 14, 1304 where The state transition matrix of the one tooth passing period is expressed as According to the Floquet theory, if any module of the eigenvalues of the state transition matrix exceeds one, the system motion is unstable; otherwise, if all modules of the eigenvalues of the state transition matrix are less than one, the system motion is stable [32]. Therefore, the boundary curve for the stable and unstable regions in the lobe diagram can be used as a criterion to judge whether chatter occurs.

Simulation Analysis
To verify the effectiveness and superiority of the proposed method in terms ofconvergence rate, prediction accuracy, and computational efficiency, a comparison between the NCM and computer numerical control milling (CNCM) methods proposed in this paper was conducted. Herein, a single-DOF benchmark example is utilized for verification, and its mathematical model [22] is written as ..
where ζ denotes the damping ratio, ω n denotes the natural frequency, and m t denotes the modal mass of the tool. The h(t) denotes the specific cutting force coefficient, which is given by Equation (10). The state-space expression of Equation (48) is the same as Equation (6), where A(t) and B(t) satisfy For an equal comparison, the parameters used in this section and listed in Table 1 are the same as those in [24], and all algorithms run on the same desktop computer (AMD Ryzen 5 5600H; CPU 4.0 GHz, 16 GB). Table 1. Physical and machining parameters of the benchmark example with a single DOF [24].

Convergence Rate
The convergence rate is regarded as an important evaluation index of algorithm accuracy. It can be expressed as the difference between the approximate and exact values. Li et al. [24] proved that the local discretization error of NCM is o h 5 . According to the discussion above, three intermediate interpolation points are considered, which makes the local discretization error of the CNCM proposed in this paper reach o h 6 . The convergence rate can be further illustrated by a single-DOF benchmark example with the following machining conditions. The different spindle speeds and depths of cut are used: Ω = 6000 rpm, a p = 0.28 mm; Ω = 6000 rpm, a p = 0.8 mm; Ω = 6000 rpm, a p = 1 mm; Ω = 10,000 rpm, a p = 0.28 mm; Ω = 10,000 rpm, a p = 0.8 mm; Ω = 10,000 rpm, a p = 1 mm. The exact value of λ0 is calculated using SDM with n = 500. Figure 1 presents the convergence rates of the NCM and CNCM proposed in this paper. The comparison results show that the CNCM proposed in this paper obtains a higher convergence rate than the NCM. In addition, the NCM suffers from more violent oscillations than the CNCM proposed in this paper. These phenomena mentioned above can illustrate that CNCM can quickly and stably reach the target value. To further analyzethe superiority of CNCM, prediction accuracy is carried out.

Prediction Accuracy
The stability lobe diagrams are given in Table 2 to compare the prediction accuracies of the NCM and CNCM proposed in this paper. The time intervals of the NCM and CNCM proposed in this paper are set to 26, 40, and 56. The reference stability boundaries denoted by the black dashed line are calculated by the SDM with n = 500. The stability lobe diagrams are constructed over a 200 × 100 grid, and the machining parameters are set as follows: Ω∈[5000 rpm, 10,000 rpm] and a p ∈ [0 mm, 5 mm]. It can be seen that the prediction accuracy of the stability lobe diagrams calculated by the CNCM proposed in this paper is higher than that calculated by the NCM. of the NCM and CNCM proposed in this paper. The comparison results show that the CNCM proposed in this paper obtains a higher convergence rate than the NCM. In addition, the NCM suffers from more violent oscillations than the CNCM proposed in this paper. These phenomena mentioned above can illustrate that CNCM can quickly and stably reach the target value. To further analyzethe superiority of CNCM, prediction accuracy is carried out.

Computational Efficiency
A stability analysis method with a high convergence rate typically requires more computation time [34]. To verify that the proposed method has a relatively ideal computation efficiency, the computational times of the NCM and CNCM at the same AMRE and MSE are discussed. As shown in Figure 3, when AMRE = 1.75%, n = 43 and t = 43.31 s are determined by the NCM, whereas n = 30 and t = 37.10 s are determined by the CNCM proposed in this paper. It is noted that the computational time can be significantly reduced with the CNCM proposed in this paper. It can also be seen that, although the two discrete numbers are different, the prediction stability boundary obtained with the CNCM proposed in this paper is below the reference stability boundary when Ω∈ [5000 rpm, 6000 rpm]. To avoid chatter, the machining parameters selected according to the above-mentioned prediction stability boundary are more reliable. In addition, when MSE = 1.12 ×

Computational Efficiency
A stability analysis method with a high convergence rate typically requires more computation time [34]. To verify that the proposed method has a relatively ideal computation efficiency, the computational times of the NCM and CNCM at the same AMRE and MSE are discussed. As shown in Figure 3, when AMRE = 1.75%, n = 43 and t = 43.31 s are determined by the NCM, whereas n = 30 and t = 37.10 s are determined by the CNCM proposed in this paper. It is noted that the computational time can be significantly reduced with the CNCM proposed in this paper. It can also be seen that, although the two discrete numbers are different, the prediction stability boundary obtained with the CNCM proposed in this paper is below the reference stability boundary when Ω∈ [5000 rpm, 6000 rpm]. To avoid chatter, the machining parameters selected according to the above-mentioned prediction stability boundary are more reliable. In addition, when MSE = 1.12 ×

Computational Efficiency
A stability analysis method with a high convergence rate typically requires more computation time [34]. To verify that the proposed method has a relatively ideal computation efficiency, the computational times of the NCM and CNCM at the same AMRE and MSE are discussed. As shown in Figure 3, when AMRE = 1.75%, n = 43 and t = 43.31 s are determined by the NCM, whereas n = 30 and t = 37.10 s are determined by the CNCM proposed in this paper. It is noted that the computational time can be significantly reduced with the CNCM proposed in this paper. It can also be seen that, although the two discrete numbers are different, the prediction stability boundary obtained with the CNCM proposed in this paper is below the reference stability boundary when Ω∈ [5000 rpm, 6000 rpm]. To avoid chatter, the machining parameters selected according to the above-mentioned prediction stability boundary are more reliable. In addition, when MSE = 1.12 ×

Computational Efficiency
A stability analysis method with a high convergence rate typically requires more computation time [34]. To verify that the proposed method has a relatively ideal computation efficiency, the computational times of the NCM and CNCM at the same AMRE and MSE are discussed. As shown in Figure 3, when AMRE = 1.75%, n = 43 and t = 43.31 s are determined by the NCM, whereas n = 30 and t = 37.10 s are determined by the CNCM proposed in this paper. It is noted that the computational time can be significantly reduced with the CNCM proposed in this paper. It can also be seen that, although the two discrete numbers are different, the prediction stability boundary obtained with the CNCM proposed in this paper is below the reference stability boundary when Ω∈ [5000 rpm, 6000 rpm]. To avoid chatter, the machining parameters selected according to the above-mentioned prediction stability boundary are more reliable. In addition, when MSE = 1.12 ×

Computational Efficiency
A stability analysis method with a high convergence rate typically requires more computation time [34]. To verify that the proposed method has a relatively ideal computation efficiency, the computational times of the NCM and CNCM at the same AMRE and MSE are discussed. As shown in Figure 3, when AMRE = 1.75%, n = 43 and t = 43.31 s are determined by the NCM, whereas n = 30 and t = 37.10 s are determined by the CNCM proposed in this paper. It is noted that the computational time can be significantly reduced with the CNCM proposed in this paper. It can also be seen that, although the two discrete numbers are different, the prediction stability boundary obtained with the CNCM proposed in this paper is below the reference stability boundary when Ω∈ [5000 rpm, 6000 rpm]. To avoid chatter, the machining parameters selected according to the above-mentioned prediction stability boundary are more reliable. In addition, when MSE = 1.12 ×

Computational Efficiency
A stability analysis method with a high convergence rate typically requires more computation time [34]. To verify that the proposed method has a relatively ideal computation efficiency, the computational times of the NCM and CNCM at the same AMRE and MSE are discussed. As shown in Figure 3, when AMRE = 1.75%, n = 43 and t = 43.31 s are determined by the NCM, whereas n = 30 and t = 37.10 s are determined by the CNCM proposed in this paper. It is noted that the computational time can be significantly reduced with the CNCM proposed in this paper. It can also be seen that, although the two discrete numbers are different, the prediction stability boundary obtained with the CNCM proposed in this paper is below the reference stability boundary when Ω∈ [5000 rpm, 6000 rpm]. To avoid chatter, the machining parameters selected according to the above-mentioned prediction stability boundary are more reliable. In addition, when MSE = 1.12 × To further verify the prediction accuracy, the arithmetic mean of the relative error (AMRE) and the mean squared error (MSE) [33] are introduced. They are defined as where a i and a i0 denote the predicted and reference axial cutting depths, respectively, and n e denotes the number of discrete points. The AMREs and MSEs are depicted in Figure 2.
The AMRE results are given below: When the time interval is 26, the AMREs of the NCM and CNCM proposed in this paper are 21.31% and 6.69%, respectively. When the time interval is 40, the AMREs of the NCM and CNCM proposed in this paper are 3.88% and 0.41%, respectively. When the time interval is 56, the AMREs of the NCM and CNCM proposed in this paper are 0.84% and 0.41%, respectively. The MSE results are given below: When the time interval is 26, the MSEs of the NCM and CNCM proposed in this paper are 2.46 × 10 −7 and 9.60 × 10 −8 , respectively. When the time interval is 40, the MSEs of the NCM and CNCM proposed in this paper are 9.85 × 10 −9 and 2.86 × 10 −10 , respectively. When the time interval is 56, the MSEs of the NCM and CNCM proposed in this paper are 1.73 × 10 −10 and 5.28 × 10 −11 , respectively. The results above demonstrate that the proposed CNCM can achieve higher prediction accuracy than the NCM.
where ai and ai0 denote the predicted and reference axial cutting depths, respectively, and ne denotes the number of discrete points. The AMREs and MSEs are depicted in Figure 2.
The AMRE results are given below: When the time interval is 26, the AMREs of the NCM and CNCM proposed in this paper are 21.31% and 6.69%, respectively. When the time interval is 40, the AMREs of the NCM and CNCM proposed in this paper are 3.88% and 0.41%, respectively. When the time interval is 56, the AMREs of the NCM and CNCM proposed in this paper are 0.84% and 0.41%, respectively. The MSE results are given below: When the time interval is 26, the MSEs of the NCM and CNCM proposed in this paper are 2.46 × 10 −7 and 9.60 × 10 −8 , respectively. When the time interval is 40, the MSEs of the NCM and CNCM proposed in this paper are 9.85 × 10 −9 and 2.86 × 10 −10 , respectively. When the time interval is 56, the MSEs of the NCM and CNCM proposed in this paper are 1.73 × 10 −10 and 5.28 × 10 −11 , respectively. The results above demonstrate that the proposed CNCM can achieve higher prediction accuracy than the NCM.

Computational Efficiency
A stability analysis method with a high convergence rate typically requires more computation time [34]. To verify that the proposed method has a relatively ideal computation efficiency, the computational times of the NCM and CNCM at the same AMRE and MSE are discussed. As shown in Figure 3, when AMRE = 1.75%, n = 43 and t = 43.31 s are determined by the NCM, whereas n = 30 and t = 37.10 s are determined by the CNCM proposed in this paper. It is noted that the computational time can be significantly reduced with the CNCM proposed in this paper. It can also be seen that, although the two discrete numbers are different, the prediction stability boundary obtained with the CNCM proposed in this paper is below the reference stability boundary when Ω∈ [5000 rpm, 6000 rpm]. To avoid chatter, the machining parameters selected according to the above-mentioned prediction stability boundary are more reliable. In addition, when MSE = 1.12 × 10 −10 , the discrete numbers and the computational time, as well as the prediction stability boundaries, are shown in Figure 4. The discrete number obtained with the NCM is greater than that obtained with the CNCM proposed in this paper. However, the computational efficiency of the NCM is slightly lower than that of the CNCM proposed in this paper. 10 −10 , the discrete numbers and the computational time, as well as the prediction stability boundaries, are shown in Figure 4. The discrete number obtained with the NCM is greater than that obtained with the CNCM proposed in this paper. However, the computational efficiency of the NCM is slightly lower than that of the CNCM proposed in this paper.

Experimental Verification
To further verify the effectiveness and superiority of the proposed CNCM, milling experiments were carried out with a feed per tooth of 0.06 mm, Ω∈ [3000 rpm, 8000 rpm] and ap∈ [0 mm, 4 mm]. The milling system was a benchmark system with two DOFs. The workpiece material was aluminum 7075. Before milling experiments, the frequency response function (FRF) of tool should be obtained first. Thus, the hammer experiment is carried out. The measuring instruments include the Kistler 9724A5000 force hammer, Endvco 2256B-10 acceleration sensor, ATI omega-160 force sensor, DEWEsoft SIRIUSI-8xACC data acquisition system and DEWEsoft data analysis system, as shown in Figure  5. The physical and machining parameters of the Harmer experiment are listed in Table 3. A stability lobe diagram of the proposed CNCM is depicted in Figure 6. It can be noted that 20 points are selected from the stability lobe diagram as the test points. 10 −10 , the discrete numbers and the computational time, as well as the prediction stability boundaries, are shown in Figure 4. The discrete number obtained with the NCM is greater than that obtained with the CNCM proposed in this paper. However, the computational efficiency of the NCM is slightly lower than that of the CNCM proposed in this paper.

Experimental Verification
To further verify the effectiveness and superiority of the proposed CNCM, milling experiments were carried out with a feed per tooth of 0.06 mm, Ω∈ [3000 rpm, 8000 rpm] and ap∈ [0 mm, 4 mm]. The milling system was a benchmark system with two DOFs. The workpiece material was aluminum 7075. Before milling experiments, the frequency response function (FRF) of tool should be obtained first. Thus, the hammer experiment is carried out. The measuring instruments include the Kistler 9724A5000 force hammer, Endvco 2256B-10 acceleration sensor, ATI omega-160 force sensor, DEWEsoft SIRIUSI-8xACC data acquisition system and DEWEsoft data analysis system, as shown in Figure  5. The physical and machining parameters of the Harmer experiment are listed in Table 3. A stability lobe diagram of the proposed CNCM is depicted in Figure 6. It can be noted that 20 points are selected from the stability lobe diagram as the test points.
The amplitude spectrum of the milling force signal can show many details in the frequency domain, which can be used to verify the prediction accuracy of the SLD presented in Figure 6. Ten critical points are selected as the milling parameters, and the milling force   Figure 5. Experiments platform, measuring instruments and calibration process. The amplitude spectrum of the milling force signal can show many deta quency domain, which can be used to verify the prediction accuracy of the SL The amplitude spectrum of the milling force signal can show many details in the frequency domain, which can be used to verify the prediction accuracy of the SLD presented in Figure 6. Ten critical points are selected as the milling parameters, and the milling force is measured by an ATI force sensor with a sample rate of 8000. When the milling process is stable, there is only one tooth passing frequency in the spectrum. However, when chatter occurs, there will be a series of disorders in the frequency spectrum. As shown in Figure 7, the milling system has the tooth passing frequency and its harmonic frequency (red value) during stable cutting situations corresponding to A1 (3500 rpm, 0.5 mm), B2 (4500 rpm, 1 mm), C1 (5500 rpm, 0.5 mm), D3 (6500 rpm, 1.5 mm), D4 (6500 rpm, 2 mm), and E1 (7500 rpm, 0.5 mm). It also has chatter frequency (blue value) during unstable cutting situations, corresponding to A2 (3500 rpm, 1 mm), B3 (4500 rpm, 1.5 mm), C2 (5500 rpm, 1 mm), and E2 (7500 rpm, 1 mm). It can be concluded that select milling parameters above the boundary will chatter, while those below the boundary can stabilize the milling process. is measured by an ATI force sensor with a sample rate of 8000. When the milling process is stable, there is only one tooth passing frequency in the spectrum. However, when chatter occurs, there will be a series of disorders in the frequency spectrum. As shown in Figure 7, the milling system has the tooth passing frequency and its harmonic frequency (red value) during stable cutting situations corresponding to A1 (3500 rpm, 0.5 mm), B2 (4500 rpm, 1 mm), C1 (5500 rpm, 0.5 mm), D3 (6500 rpm, 1.5 mm), D4 (6500 rpm, 2 mm), and E1 (7500 rpm, 0.5 mm). It also has chatter frequency (blue value) during unstable cutting situations, corresponding to A2 (3500 rpm, 1 mm), B3 (4500 rpm, 1.5 mm), C2 (5500 rpm, 1 mm), and E2 (7500 rpm, 1 mm). It can be concluded that select milling parameters above the boundary will chatter, while those below the boundary can stabilize the milling process.  Figure 8. When chatter occurs, the quality of the machined surface is poor. It can be seen that the chatter becomes increasingly obvious as the depth of the cut increases. The machined surface roughness values are listed in Table 4. Under the conditions of A1-A4, the average values of the machined surface roughness are 1.358, 2.392, 2.462, and 2.910 µm, respectively. It is obvious that the chatter significantly decreases the machined surface quality. The machined surface roughness under the conditions of B1-B4, C1-C4, D1-D4, and E1-E4 is also listed in Table 4. It can be observed that, under the conditions of A1, B1, B2, C1, D1-D4, and E1, the average values of the machined surface roughness are 1.358, 1.181, 1.501, 1.363, 1.312, 1.514, 1.655, 1.937, and 1.305 µm, respectively. The results show that the cutting parameters below the stability boundary can avoid chatter as much as possible, thereby ensuring the quality of the machined surface. The machined surfaces corresponding to A1-A4 are shown in Figure 8. When chatter occurs, the quality of the machined surface is poor. It can be seen that the chatter becomes increasingly obvious as the depth of the cut increases. The machined surface roughness values are listed in Table 4. Under the conditions of A1-A4, the average values of the machined surface roughness are 1.358, 2.392, 2.462, and 2.910 µm, respectively. It is obvious that the chatter significantly decreases the machined surface quality. The machined surface roughness under the conditions of B1-B4, C1-C4, D1-D4, and E1-E4 is also listed in Table 4

Conclusions
This study focuses on the convergence rate and prediction accuracy of the milling stability based on the composited Newton-Cotes formula. For example, when the discrete number is 40, the convergence rate and prediction accuracy improve about 4.1 and 9.4 times, respectively, while computational efficiency is similar. The reason is that because the proposed method considers three intermediate interpolation points, the local discrete error of the proposed method is smaller than that of the comparison method, which ensures the prediction accuracy of the proposed method. The state transition matrix of the proposed method is divided into three submatrices, which can avoid double calculation and matrix multiplication. The simulation results of the milling model with a single DOF validate that the proposed method has a high convergence rate and prediction accuracy while ensuring computational efficiency. The experimental results of the milling model with two DOFs validate that the machining parameters below the prediction stability boundary obtained with the proposed method can avoid chatter as much as possible while ensuring the machined surface quality.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.