## 1 Introduction  
The theory of variations has been used so far to derive integrators that preserve the physical properties of the system in question, e.g., a symplectic form [1, 2, 3, 4]. When tested on several Hamiltonian ordinary differential equations (ODEs), symplectic integrators have shown robustness, efficiency and accuracy even for long integrations [5, 6, 7, 8, 9, 10]. Hamiltonian systems that lead to partial differential equations (PDEs) allow a similar description, referred to as the \"Veselov discretisation for PDEs in variational form\" [3, 4]. Due to the presence of more preserved quantities, and hence more symplectic forms, these integrators (named multisymplectic schemes) have been derived on the basis of the structure of Hamiltonian PDEs from Lagrangian or Hamiltonian perspectives [11, 12, 13, 14].  
In our recent works [15, 16] we investigated the weaknesses (e.g., numerical instabilities) of standard multisymplectic integrators derived via simple quadratic expressions adopted for the spatial and time derivatives [11]. We also explored several schemes that use complex expressions based on the non-standard finite differences of [17]. The latter proved to be more accurate, and their performances are better than those the standard ones. They use trigonometric functions in the derivative expressions and preserve certain properties and structures dictated by the physical system [18, 19, 16].  
In this work, we extend our previous multisymplectic schemes by exploiting the advantages of exponential functions to express the spatial and time derivatives. The focus is on the solution of the wave equation, a hyperbolic-type equation within the theory of linear PDEs. Unlike other PDEs where the solutions are smooth when the initial conditions are smooth, the hyperbolic equations are challenging due to their sharp solutions behaviour [20, 21].  
In addition, we consider a non-linear version of the standard wave equation, namely, the non-linear seismic wave equation [22]. To that end, we obtain its multisymplectic form directly from the variational principle, i.e., by using the Lagrangian approach. For simplicity, in this article we work in local coordinates, although a more general coordinate system can be considered. Our attempt leads to multi-symplectic integrators with very good energy performance with respect to conservation of the nearby Hamiltonian, up to an exponentially small error [11; 15; 16; 23; 24; 25; 26].  
The paper is organised as follows. Firstly, we introduce some notions from the smooth multisymplectic Hamiltonian systems in Section 2. These are used in Section 3 to formulate a discrete version of the Lagrangian density relying on exponential functions and triangular space-time discretisation (a square space-time discretisation is additionally presented in the Appendix A). The resulting discrete Lagrangian is then applied in the discrete version of the Euler-Lagrange field equations in order to deduce the corresponding discrete equations. The proposed scheme is tested in the linear wave equation in Section 4 and in the non-linear seismic wave equation in Section 5. Subsequently, error analysis and convergence plots are presented to demonstrate the improved numerical properties and the behaviour of the proposed scheme as compared to the standard ones. Finally, the main conclusions extracted from our study are summarised in Section 6.