Summary: In this paper we develop and analyze a new superconvergent local discontinuous Galerkin (LDG) method for approximating solutions to the fourth-order Euler-Bernoulli beam equation in one space dimension. We prove the \(L^2\) stability of the scheme and several optimal \(L^2\) error estimates for the solution and for the three auxiliary variables that approximate derivatives of different orders. Our numerical experiments demonstrate optimal rates of convergence. We also prove superconvergence results towards particular projections of the exact solutions. More precisely, we prove that the LDG solution and its spatial derivatives (up to third order) are \(\mathcal O (h^{k+3/2})\) super close to particular projections of the exact solutions for \(k\)th-degree polynomial spaces while computational results show higher \(\mathcal O (h^{k+2})\) convergence rate. Our proofs are valid for arbitrary regular meshes and for \(P^k\) polynomials with \(k\geq 1\), and for periodic, Dirichlet, and mixed boundary conditions. These superconvergence results will be used to construct asymptotically exact a posteriori error estimates by solving a local steady problem on each element. This will be reported in Part II of this work, where we will prove that the a posteriori LDG error estimates for the solution and its derivatives converge to the true errors in the \(L^2\)-norm under mesh refinement.