The novel finite element formulations fall into the category of geometrically exact Kirchhoff-Love beams. A prominent characteristic of this category is that the absence of shear deformation is strongly enforced by removing two degrees of freedom. Further, the corresponding beam theories exhibit not only translational but also rotational degrees of freedom and their configurations thus form a non-additive and non-commutative space. Sophisticated interpolation schemes are required that need to be tested not only for locking, spatial convergence behavior, and energy conservation, but also for observer invariance and path-independence. For the three novel beam element formulations all these properties are analytically and numerically studied and confirmed, if applicable. Two different rotation parameterization strategies are employed based on the well-known geodesic interpolation used in many Simo-Reissner beams and the lesser known split into the so-called \textit{smallest rotation} and a torsional part. Application of the former parameterization results in a mixed finite element formulation intrinsically free of locking phenomena. Additionally, the first geometrically exact Kirchhoff-Love beam element is presented, which strongly enforces inextensibility by removing another degree of freedom. Furthermore, the numerical efficiency of the new beam formulations is compared to other beam elements that allow for or suppress shear deformation. When modeling very slender beams, the new elements offer distinct numerical advantages. Standard molecular dynamics simulations, which are commonly used to study polymers, suffer from a lack of a careful mathematical basis and the use of an expensive explicit time integration scheme. To circumvent these shortcomings and to be able to simulate stretching experiments on relevant time scales, the problem is described by a stochastic partial differential equation, which can be solved using the finite element method with a backward Euler temporal discretization. In detail, the polymer is represented by a Kirchhoff-Love beam with a linear elastic constitutive model. Inertial and electrostatic forces are neglected. It is deformed by a distributed load mimicking collisions with molecules of the surrounding fluid. Naturally, this load heavily fluctuates over time and space and mean values need to be computed in a Monte Carlo manner. To vastly speed up the fitting process to experimental data in a Bayesian framework, a surrogate model based on a Gaussian process is set up, which directly computes the mean values for given material parameters. The uncertainties and correlations of the material parameters are studied and compared to the literature.