A spatial discretization and substructure method is developed to accurately calculate dynamic responses of one-dimensional structural systems, which consist of length-variant distributed-parameter components, such as strings, rods, and beams, and lumped-parameter components, such as point masses and rigid bodies. The dependent variable of a distributed-parameter component is decomposed into boundary-induced terms and internal terms. The boundary-induced terms are interpolated from boundary motions, and the internal terms are approximated by an expansion of trial functions that satisfy the corresponding homogeneous boundary conditions. All the matching conditions at the interfaces of the components are satisfied, and the expansions of the dependent variables of the distributed-parameter components absolutely and uniformly converge if the dependent variables are smooth enough. Spatial derivatives of the dependent variables, which are related to internal forces/moments of the distributed-parameter components, such as axial forces, bending moments, and shear forces, can be accurately calculated. Combining component equations that are derived from Lagrange's equations and geometric matching conditions that arise from continuity relations leads to a system of differential algebraic equations (DAEs). When the geometric matching conditions are linear, the DAEs can be transformed to a system of ordinary differential equations (ODEs), which can be solved by an ODE solver. The methodology is applied to several moving elevator cable-car systems in Part II of this work.