## 1 Introduction
### Spectral System Equation and Dispersion Relation
If there is not, directly go to step (6); If there is, search the root by virtue of the dichotomizing method through inspecting the real or imaginary parts of the determinants or finish the root searching by the golden section method by observing the amplitudes of the determinants. * Judge if adequate numbers of the requested quantity have been obtained or judge if the second tentative value is exceeding the required range. If so, go to step (7). If not, go to the next step. * Assign the second tentative value as the first one and then return to step (2). * Judge if adequate numbers of root have been obtained for the current value of known quantity or judge if the current value of known quantity exceeds the required range. If not, start the calculation from step (1) again. If yes, terminate the calculation. After all the roots within the required range of the requested quantity have been obtained for all the specified values of known quantity with a given range, the dispersion spectra between the known and requested quantities in the corresponding area can be plotted. These dispersion curves represent the dispersion properties of the characteristic waves from a viewpoint. It should be emphasized that the involved exponential functions in system Equation (15) and the dispersion Equation (16) are \(\mathrm{e}^{\mathrm{i}ql}\) from the global scattering relation (12) and \(\mathrm{e}^{-\mathrm{i}k_{r(j)}l_{(j)}}\) from the global phase relation (A3), respectively. Since the elaborately chosen \(k_{r(j)}\) has guaranteed \(\mathrm{Re}[-\mathrm{i}k_{r(j)}l_{(j)}]\leq 0\), thus \(\mathrm{e}^{-\mathrm{i}k_{r(j)}l_{(j)}}\) is always an exponentially decaying function. Comparatively, the other analytical methods like the well-known transfer matrix method (TMM) contain both \(\mathrm{e}^{-\mathrm{i}k_{r(j)}l_{(j)}}\) and \(\mathrm{e}^{\mathrm{i}k_{r(j)}l_{(j)}}\) in its formulation, and thus inevitably both exponentially decaying and growing functions are involved, which will lead to a numerically ill-condition problem as discussed in Ref. [79]. Our MRRM formulation, however, should be numerically well-conditioned as long as \(\mathrm{Re}[\mathrm{i}ql]\leq 0\) satisfies or \(\mathrm{Re}[\mathrm{i}ql]\) is not too big. These conditions are relatively easy to guarantee since in the wave problem studied here the eigenvalue \([\mathrm{i}ql]\) is symmetrical with respect to \(q=0\), i.e., the characteristic waves in opposite directions should give the same dispersion property. Thus, our proposed MRRM is applicable for accurate calculation of the dispersion curves until a relatively high frequency where the Timoshenko beam theory still holds for modeling structural members. With regard to the aforementioned numerically well-conditioned property, the advantage of our MRRM as compared to TMM will be illustrated by numerical examples in Section 3.2. It should also be pointed out that our proposed MRRM formulation based on the Timoshenko beam theory in this section can actually be degenerated to a counterpart formulation for the bi-coupled periodic Euler-Bernoulli beam structures, as long as \(c_{t}=\infty\), \(\beta_{r}=\mathrm{i}k_{r}\), \(\zeta_{r}=-(\mathrm{i}k_{r})^{3}EI\), \(\mu_{r}=(\mathrm{i}k_{r})^{2}EI\), \(k_{2}=k=\sqrt[4]{\omega^{2}A/(c_{t}^{2}I)}\) and \(k_{3}=-\mathrm{i}k_{2}=-\mathrm{i}k\) are adopted [76] in all the above related formulas. Again, the MRRM formulation in this case is valid as long as the Euler-Bernoulli beam theory is effective.